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EXECUTIVE  SUMMARY 


OBJECTIVE 


The  purpose  of  this  study  was  to  demonstrate  the  feasibility  of  using  a  small  computer  * 
to  operate  a  three-dimensional  atmospheric  model  to  forecast  the  transport  and  dispersion  of 
airborne  hazardous  materials  at  Vandenberg  Air  Force  Base  (VAFB), 

BACKGROUND 

Air  Force  (AF)  models,  designed  to  calculate  atmospheric  transport  and  dispersion  at  VAFB, 
assume  that  the  wind  field  over  the  geographical  area  of  interest  can  be  adequately  represented 
by  the  wind  speed  and  direction  measured  at  a  single  point  or  by  a  vertical  profile  of  wind 
speed  and  direction  over  a  single  site.  While  such  an  assumption  is  acceptable  over  flat  terrain 
and  over  a  short  time,  the  complex  terrain  and  meteorology  at  VAFB  and  the  reduction  in 
personal  exposure  limits  for  the  AF  chemicals  of  interest  require  that  more  accurate  methods 
be  used  to  make  these  calculations. 

Current  efforts  to  improve  VAFB’s  modeling  capability  also  assume  the  three-dimensional 
wind  field  can  be  represented  by  a  single  wind  value.  The  complexity  of  the  terrain  and 
meteorology  at  VAFB  require  a  model  that  takes  into  account  the  three-dimensional  wind  field 
over  the  entire  VAFB  region. 

Recent  advancements  in  computer  hardware  technology,  however,  have  the  potential  to 
allow  much  more  sophisticated  models  to  be  operated  in  much  less  time  than  previously 
considered  possible.  The  contractor  investigated  the  feasibility  of  using  the  three-dimension^ 
atmospheric  models  HOTMAC  (Higher  Order  Turbulence  Model  for  Atmospheric  Circulations) 
and  RAPTAD  (Random  Puff  Transport  and  Diffusion)  to  meet  the  objective  of  the  project. 

TASKS 

The  specific  tasks  are  described  in  Section  IL  The  first  task  was  to  prepare  input  files  for 
HOTMAC  and  RAPTAD.  The  input  data  included  digitized  terrain  elevations,  initial  wind  and 
potential  temperature  distributions,  and  the  surface  condition.  The  main  task  was  to  operate 
HOTMAdVRAPTAD  with  VAFB  terrain  data  and  compare  output  with  the  observations  and 
that  obtained  by  AF  models.  The  contractor  operated  the  HOTMAC/RAPTAD  model  on 
a  Micro Vax  and  a  Sun  work  station.  Finally,  the  contractor  investigated  the  feasibility  of 
incorporating  positive  and  negative  source  buoyancy  effects  in  HOTMAC  and  RAPTAD. 

METHODOLOGY 

Section  in  reviewed  briefly  the  model  equations,  boimdary  conditions,  and  numerical 
schemes  used  in  HOTMAC  and  RAPTAD.  A  detailed  model  description  is  given  in  Appendix  A. 

HOTMAC  is  a  mesoscale  atmospheric  model  that  can  forecast  three-dimensional  distri¬ 
butions  of  wind  speed,  wind  direction,  turbulence,  temperature,  and  water  vapor.  The  basic 


equations  for  HOTMAC  are  the  conservation  equations  for  mass,  momentum,  internal  energy, 
mixing  ratio  of  water  vapor,  and  turbulence  kinetic  energy. 

RAPTAD  is  a  Lagrangian  puff  code  based  on  the  Monte  Carlo  statistical  diffusion  process. 
The  center  location  and  standard  deviation  of  concentration  distribution  for  each  puff  are 
computed  by  use  of  wind  and  turbulence  modeled  by  HOTMAC.  Then,  concentration  at  any 
location  is  computed  by  summation  of  concentrations  contributed  by  all  the  puffs.  RAPTAD  can 
be  used  under  extreme  conditions  with  highly  heterogeneous  wind  and  turbulence  distributions 
where  a  convectional  Gaussian  plume  model  may  fail. 

TEST  DESCRIPTION 

To  test  the  feasibility  of  operating  HOTMAC/RAPTAD  with  the  VAFB  terrain,  we  seiec:..  J 
data  from  the  Mountain  Iron  (MI)  diffusion  experiments,  which  were  conducted  at  VAFB 
during  1965  and  1966.  The  objective  of  the  MI  program  was  to  establish  quantitative  diffusion 
predictions  for  use  as  range  safety  tools  in  the  South  Vandenberg  (SV)  ballistic  and  space 
vehicle  operations  area. 

Meteorological  masts  were  deployed  to  measure  wind  speeds  and  wind  directions.  These 
wind  data  were  used  to  examine  the  performance  of  HOTMAC.  Vertical  profiles  of  wind  speed, 
wind  direction,  and  wet  and  dry  temperature  up  to  2000  meters  above  the  ground  were  measured 
by  rawinsondes  at  four  sites,  and  wind  profiles  were  taken  by  a  wiresonde  at  the  tracer  release 
site.  The  vertical  profiles  of  winds  and  temperatures  were  used  to  initialize  the  model  variables 
and  test  the  model  predictions. 

Fluorescent  particles  were  released  to  aid  the  understanding  of  transport  and  diffusion 
processes  and  the  derivation  of  an  empirical  formula  for  the  pollutant  concentration  distribution 
in  the  SV  area. 

A  total  of  113  tracer  releases  were  made,  from  which  we  selected  the  MI87,  MI90,  and 
MI91  cases  to  evaluate  the  performance  of  HOTMAC  and  RAPTAD. 

RESULTS 

MI87  Simulation 

The  meteorological  condition  for  the  MI87  case  was  characterized  as  “weak  ambient  and 
sea-breeze  flow,  with  little  distinction  between  wind  vectors,  but  sharp  cooling  with  fog  rolling 


The  initial  HOTMAC  potential  temperature  profile  was  determined  by  averaging  four 
upper  air  soundings  at  1300  local  standi  time  Ost).  The  potential  temperature  profile  was 
approximately  0.044°  C/m  from  the  sea  surface  to  460  meters  above  mean  sea  level  (msl), 
0.0142  °C/m  between  460  meters  and  960  meters  above  msl,  and  0.0045°C/m  above  960 
meters  above  msl.  Initial  wind  speed  and  direction  were  determined  by  examining  five  upper 
air  soundings  at  1300  1st.  Initial  upper  air  wind  speed  and  wind  direction  above  the  boundary 
layer  were  estimated  to  be  2  m/s  and  225  degrees,  respectively. 


IV 


The  computational  domain  was  40  x  48  km^  with  a  horizontal  grid  spacing  of  1  km.  To 
resolve  the  details  of  topography  in  the  vicinity  of  the  tracer  release  site,  we  nested  a  fine- 
resolution  grid  15  X  16  km^  with  a  horizontal  grid  spacing  of  0.5  km. 

Integration  started  at  0500  1st,  June  13,  1966,  and  continued  for  over  12  hoiu^.  The  plume 
was  released  at  1310 1st  for  30  minutes  as  was  done  in  the  experiment  The  plume  was  followed 
for  4  hours  in  the  model  computation.  By  that  time,  the  plume  was  transported  far  away  firom 
the  sampling  areas. 

The  modeled  horizontal  wind  vectors  at  10  meters  above  the  ground  at  1300  1st,  June  13, 
1966  were  in  good  agreement  with  the  observed  winds  in  the  surface  layer,  although  observed 
wind  vectors  show  considerable  variation  in  direction  and  magnitude  compared  with  modeled 
winds. 

Less  satisfactory  results  were  obtained  when  vertical  profiles  of  the  modeled  wind  speed 
and  wind  direction  were  compared  with  observations.  Wind  speed  and  wind  direction  become 
highly  variable  in  space  and  time  when  wind  speed  is  small.  It  should  be  noted  that  the 
observations  were  instantaneous  values  whereas  the  modeled  results  were  ensemble  averages. 

The  ground-level  exposure  values  were  computed  by  integrating  concentration  with  time 
until  the  concentration  became  negligible.  Although  tiie  modeled  plume  direction  did  not 
exactly  match  the  observed  one,  the  modeled  ground-level  exposures  along  the  plume  axis 
were  in  good  agreement  with  observations. 


MI90  Simulation 

Meteorological  conditions  for  the  MI90  case  were  significantly  different  firom  those  for  the 
MI87  case.  The  prevailing  wind  speed  for  MM)  was  close  to  13.5  m/s  compared  with  2  m/s  for 
MI87.  In  addition,  for  MI90  a  deep,  well-mixed  layer  existed  for  the  first  500  meters  above  the 
ground,  which  probably  resulted  from  strong  mixing  caused  by  turbulence  generated  by  high 
winds.  The  potential  temperature  profiles  for  MI87  exhibited  shallow  (100  meters)  mixed-layer 
depths  that  are  typical  on  the  west  coast  as  a  result  of  large-scale  subsidence. 

The  initial  potential  temperature  profile  was  determined  from  the  sounding  at  the  Building 
22  site.  The  potential  temperature  lapse  rate  was  almost  neutral  (0.0002  ®C/m)  from  the  surface 
to  670  meters  above  msl,  0.0611  ®C/m  between  670  meters  and  880  meters  above  msl,  and 
0.0070  ®C/m  above  880  meters  above  msl.  Initial  wind  speed  (13.5  m/s)  and  direction  (332 
degrees)  were  determined  by  examination  of  three  upper  air  soundings. 

Integration  started  at  0600  1st,  June  21,  1966  and  continued  for  24  hours  until  0600  1st  the 
following  day.  The  plume  was  released  for  30  minutes  starting  at  2300 1st,  June  21,  1966.  The 
plume  was  sampled  for  2  hours  in  the  model  computation.  By  this  time  the  plume  had  left  the 
computational  domain,  and  the  surface  concentration  values  were  practically  zero. 

With  the  wind  and  turbulence  distributions  computed  by  HOTMAC,  RAPTAD  simu¬ 
lated  the  groimd-level  exposure  values  which  were  in  good  agreement  with  the  observa¬ 
tions.  The  observed  exposures  at  1870  and  5620  meters  fi-om  the  source  were  5.08  x  10~^ 


and  4.97  x  10  ^  s/m^,  respectively.  The  corresponding  modeled  values  at  1895  and  5603 
meters  from  the  source  were  1.83  x  lO”^  and  3.45  x  10"^  s/m^,  respectively. 


MI91  Simulation 

The  MI91  release  was  at  0203  1st,  June  22,  1966  which  was  only  3  hours  after  the  MI90 
release.  However,  the  wind  direction  changed  considerably  to  northerly,  as  much  as  30  degrees 
from  the  MI90  case.  Based  on  the  MI91  data  log,  the  wind  direction  veering  from  northwest 
to  north  occurred  around  midnight  June  21,  1966,  and  was  most  evident  in  the  layers  between 
300  and  1300  meters  above  msl. 

Wind  directions  at  the  grid  levels  between  300  and  1300  meters  above  msl  were  nudged 
toward  the  observed  wind  directions.  We  started  the  HOTMAC  computation  with  wind  direction 
nudging,  using  the  output  at  2300  from  the  MI  90  simulation  to  initialize  the  model  variables. 

The  modeled  horizontal  wind  vectors  clearly  indicated  that  the  wind  directions  shifted  from 
the  corresponding  values  for  the  MI90  case.  The  nudging  was  not  applied  to  the  layers  in 
the  first  300  meters  above  the  ground  but  the  surface  layer  winds  were  clearly  influenced  by 
the  nudging.  This  good  commurtication  between  the  surface  and  upper  level  winds  apparently 
resulted  from  the  strong  vertical  mixing  resulting  from  extraordinary  large  turbulence  for  a 
nocturnal  period.  The  high  prevailing  wind  speed  (13.5  m/s)  was  responsible  for  the  large 
turbulence  values. 

Model  Comparisons 

There  were  a  few  attempts  in  the  past  to  simulate  the  MI  wind  and  concentration  data. 
Himter  (Reference  12)  investigated  the  performance  of  a  diagnostic  wind  model  WOCSS  for 
eight  cases  of  the  MI  data.  Thykier — ^Nielsen  et  al.  (Reference  13)  tested  a  diagnostic  wind 
model  (LINCOM)  and  a  puff  model  (RIMPUFF)  for  two  cases  of  the  MI  data.  The  data 
common  to  both  studies  were  the  MI87  and  MI90  cases.  Therefore,  we  used  these  two  cases  to 
compare  the  performances  of  WOCSS,  LINCXJM,  RIMPUFF,  HOTMAC,  and  RAPTAD,  and 
discussed  the  strengths  and  weaknesses  of  each  model. 


Application  on  Desktop  Computers 

Recent  advances  in  desktop  computer  capabilities,  particularly  those  of  an  engineering 
workstation,  are  astonishing.  A  high-p^ormance  workstation  has  reportedly  exceeded  a  super¬ 
computer  in  certain  scalar  operations.  The  affordability  and  portability  of  the  desktop  computer 
have  opened  the  door  to  many  applications  that  were  previously  considered  impossible. 

One  logical  application  is  upgrading  toxic-hazard  modeling  capabilities  for  emergency 
response  management  Scientists  at  Los  Alamos  National  Laboratory  are  in  the  process  of 
developing  three-dimensional  forecasting  models  that  run  on  a  desktop  computer. 

In  this  section  we  discussed  our  experience  in  running  HOTMAC  and  RAPTAD  on  a  Sun 
Microsystems  workstation  and  a  MictoVax  computer. 
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Buoyant  Plumes 


In  this  section,  we  investigated  the  feasibility  of  incorporating  positive  and  negative  source 
buoyancy  effects  in  HOTMAC  and  RAPTAD. 

The  most  widely  used  method  to  treat  a  buoyant  plume  is  to  use  a  plume  rise  theory 
advanced  by  Briggs  (Reference  14).  One  can  estimate  the  height  where  the  plume  levels  off. 
This  height  is  referred  to  as  an  effective  stack  height  Therefore,  the  simplest  way  to  deal  with 
a  buoyant  plume  is  to  apply  the  Brigg’s  formula  in  the  initial  stage  when  the  buoyancy  effect 
is  significant.  When  the  plume  becomes  neutrally  buoyant  RAPTAD  can  be  applied  as  before 
except  that  a  stack  height  is  replaced  by  an  effective  stack  height. 

Our  approach  is  quite  different  from  that  of  a  plume  rise  model.  We  specify  a  heat  source  as 
a  boundary  condition,  then  the  governing  equations  in  HOTMAC  compute  wind,  temperature, 
and  turbulence  distributions.  RAPTAD  uses  these  HOTMAC  variables  to  determine  the  centroid 
location  and  size  of  each  puff.  Unlike  a  plume  rise  model,  a  puff  in  RAPTAD  is  a  tracer  that 
follows  precisely  the  turbulent  fluid  motion.  Our  approach,  in  this  smdy,  is  referred  to  as  a 
“dynamic  plume  model.”  The  advantage  of  a  dynamic  plume  model  is  that  the  fluid  motions  and 
turbulent  processes  are  consistent  with  the  physical  laws  expressed  by  the  governing  equations. 
In  other  words,  we  can  minimize  the  ambiguities  and  inconsistencies  introduced  in  a  passive 
plume  model.  For  example,  methods  used  in  a  passive  plume  model  to  parameterize  turbulent 
entrainment  processes  are  based  on  simple  and  idealized  conditions.  A  passive  plume  model, 
by  definition  carmot  provide  feedback  to  the  ambient  flow  coniitions.  For  example,  vertical 
motions  resulting  from  a  plume  rise  do  not  affect  the  distribution  of  horizontal  wind  components 
in  the  background  flow.  Of  course,  this  is  inconsistent  with  mass  continuity. 

CONCLUSIONS 

•  It  is  feasible  to  operate  the  three-dimensional  atmospheric  models  HOTMAC  and 
RAPTAD  to  forecast  the  transport  and  dispersion  of  airborne  materials  at  VAFB. 

•  HOTMAC  and  RAPTAD  predictions  were  at  least  as  good  as  those  obtained  by  diagnostic 
models  and  the  MI  model  where  wind  data  were  available. 

•  HOTMAC  and  RAPTAD  predictions  were  far  better  than  those  obtained  by  diagnostic 
models  and  they  were  the  best  practical  solution  where  wind  data  were  not  available. 

•  It  is  feasible  to  operate  HOTMAC  and  RAPTAD  on  a  desktop  computer.  HOTMAC  took 
4  hours  11  minutes  and  22  hours  10  minutes  CPU  time,  respectively,  on  a  Sun  4/110 
workstation,  and  a  MicroVax  2000  computer  for  a  28  hours  forecast  with  21  x  25  x  16 
grid  points.  RAPTAD  took  26  minutes  and  3  hours  45  minutes  CPU  time,  respectively, 
on  a  Sun  4/1 10  workstation,  and  a  MicroVax  2000  computer  for  a  20  hour  simulation. 

•  The  affordability  and  portability  of  the  desktop  computer  has  opened  the  door  to 
upgrading  toxic-hazard  modeling  capabilities  for  emergency  response  management  at 
VAFB. 
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•  The  HOTMAC  and  RAPTAD  modeling  system  would  be  a  useful  addition  to  enhance 
current  VAFB  emergency  response  management  capabilities. 

•  It  is  feasible  to  incorporate  the  positive  and  negative  source  buoyancy  effects  in  HOT¬ 
MAC  and  RAPTAD.  Our  approach  can  minimize  the  ambiguities  and  inconsistencies 
introduced  in  a  simple  plume  rise  model. 

RECOMMENDATIONS 

•  Become  familiar  with  background  theories  and  operating  procedures  of  HOTMAC  and 
RAPTAD.  Install  the  models  on  a  mainframe  computer  and  repeat  the  simulations 
discussed  in  this  report. 

•  Consider  upgrading  computer  capabilities  for  emergency  response  applications  at  VAFB. 
Micro Vax  computers  may  not  be  fast  enough  to  operate  HOTMAC  and  RAPTAD. 

•  Modify  HOTMAC  to  run  all  the  time  and  store  wind  and  turbulence  data  for  the  next  24 
hours.  When  an  emergency  occurs,  RAPTAD  can  be  used  immediately  with  the  wind 
data  stored  on  a  disk.  The  RAPTAD  computation  is  much  faster  than  that  of  HOTMAC. 
Thus,  this  approach  meets  better  the  time  constraint  imposed  under  emergency  situations. 

•  Develop  a  method  to  integrate  tower  data  into  HOTMAC.  The  wind  data  can  be  used  to 
initialize  and  correct  the  wind  distribution  in  HOTMAC.  These  may  be  accomplished  by 
a  dynamical  initialization  technique  and  a  four-dimensional  data  assimilation  method. 

•  Develop  a  method  to  predict  concentration  variances  which  can  be  used  to  estimate 
uncertainties  associated  with  predictions.  A  short  time  averaging  value  is  required 
for  predicting  concentration  of  toxic  materials.  Such  a  value  normally  exhibits  great 
variations  in  time  and  space. 

•  Add  model  physics  necessary  to  simulate  the  evolution  of  fog  formation  and  dissipation 
processes.  Fog  is  frequently  observed  at  VAFB  and  it  affects  the  heat  energy  bdance 
at  the  ground. 

•  Continue  to  investigate  the  feasibility  of  incorporating  positive  and  negative  source 
buoyancy  effects  in  HOTMAC  and  RAPTAD  and  test  the  scheme  with  observations. 
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SECTION  I 


INTRODUCTION 

Air  Force  (AF)  models  designed  to  calculate  the  atmospheric  transport  and  dispersion 
of  pollutants  at  Vandenberg  Air  Force  Base  (VAFB)  include  the  Rocket  Exhaust  Effluent 
Dispersion  Model  (REEDM)  system  (Reference  1;  Reference  2)  and  the  Ocean  Breeze/Dry 
Gulch  (OB/DG)  equation  (Reference  3),  Mountain  Iron  (MI)  equation  (Reference  4;  Reference 
5),  and  Sudden  Ranch  (SR)  equation  (Reference  6).  REEDM  was  designed  to  model  the 
transport  and  dispersion  of  buoyant  plumes  of  rocket  exhaust  during  normal  launch  and  fireball 
rise  and  subsequent  dispersion  in  the  event  of  an  explosive  accident  OB/DG,  MI,  and  SR  are 
empirical  equations  designed  to  calculate  toxic  hazard  corridors  resulting  from  accidental  spills 
of  toxic  chemicals.  These  models  assume  that  the  wind  field  over  the  geographical  area  of 
interest  can  be  adequately  represented  by  the  wind  speed  and  direction  measured  at  a  single 
point  (or,  in  the  case  of  REEDM,  a  vertical  profile  of  wind  speed  and  direction  over  a  single 
point).  These  models  also  assume  that  the  toxic  plume  or  puff  released  to  the  atmosphere  is 
uniformly  distributed  around  a  centerline  determined  from  the  wind  direction  measured  at  die 
single  observation  site. 

While  the  model  and  equations  listed  above  have  been  acceptable  in  the-past,  the  complex 
terrain  and  meteorology  at  VAFB,  increased  amounts  of  toxic  chenucals  stared  there,  increases 
in  the  surrounding  civilian  population,  higher  anticipated  latmch  rates,  and  the  reduction  in 
personal  exposure  limits  for  the  AF  chemicals  of  interest  (nitrogen  tetroxide  and  Aerozine-50) 
require  more  accurate  methods.  Current  efforts  to  improve  VAFB’s  modeling  capability  include 
development  of  the  Air  Force  Toxic  (AFTOX)  dJhemical  Dispersion  model  (Reference  7)  by  the 
Air  Force  Geophysics  Laboratory  (AFGL).  AFTOX,  in  its  current  form,  however,  also  assumes 
the  three-dimensional  wind  field  can  be  represented  by  a  single  wind  value.  Efforts  have  been 
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made  to  couple  AFTOX  with  a  surface  layer  wind  flow  model;  however,  the  complexity  of  the 
terrain  and  meteorology  at  VAFB  require  a  model  that  takes  into  account  the  three-dimensional 
wind  field  over  the  entire  VAFB  region. 

Previous  AF  efforts  to  improve  toxic  hazard  models  have  been  limited  by  available 
computer  power.  Recent  developments  in  computer  hardware  technology,  however,  have  the 
potential  to  allow  much  more  sophisticated  models  to  be  operated  in  much  less  time  than 
previously  considered  possible. 

A.  OBJECTIVE 

The  purpose  of  this  study  is  to  demonstrate  the  feasibility  of  using  a  small  computer  to 
operate  three-dimensional  hydrodynamic  and  diffusion  models  to  describe  the  transport  and 
dispersion  of  atmospheric  pollutants  at  VAFB.  This  effort  will  use  the  Los  Alamos  National 
Laboratory  (LANL)  HOTMAC  (Higher  Order  Turbulence  Model  for  Atmospheric  Circulations) 
and  RAPTAD  (Random  Puff  Transport  and  Diffusion)  models  (Reference  8)  on  a  desktop 
computer  to  accomplish  the  objective. 

B.  BACKGROUND 

HOTMAC  and  RAPTAD  are  significantly  different  from  any  AF  models  mentioned  above. 
HOTMAC  is  referred  to  as  a  prognostic  model  and  solves  a  set  of  time-dependent  physical 
equations  such  as  conservation  equations  of  momentum,  intranal  energy,  mixing  ratio  of  water 
vapor,  and  turbulence  variables.  HOTMAC  forecasts  three-dimensional  distributions  of  wind 
speed,  wind  direction,  temperature,  and  moisture.  The  current  AF  models  use  a  single  wind 
value  or  a  vertical  profile  of  wind  speed  and  direction  at  a  single  location.  Clearly,  a  single  value 
or  a  single  profile  of  wind  caimot  provide  realistic  wind  distributions  for  diffusion  computations 
because  the  t^ain  and  meteorology  in  the  VAFB  area  are  highly  heterogeneous. 
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HOTMAC  provides  to  RAPTAD  both  mean  and  turbulence  variables  to  simulate  transport 
and  diffusion  processes  of  airborne  materials.  Only  a  few  mesoscale  atmospheric  models  can 
forecast  three-dimensional  variations  of  atmospheric  turbulence.  Other  mesoscale  models  are 
being  modified  to  incorporate  second-moment  turbulence  closure  equations  similar  to  those  first 
implemented  in  HOTMAC.  HOTMAC  and  RAPTAD  predictions  have  been  tested  extensively 
with  observations  and  they  have  been  running  successfully  on  a  desktop  computer.  Therefore, 
HOTMAC  and  RAPTAD  offer  a  highly  advanced  forecasting  capability,  compared  with  simple 
emergency  response  management  models. 

An  operational,  three-dimensional,  predictive  modeling  system  that  meets  all  the  require¬ 
ments  for  the  emergency  response  management  at  VAFB  is  a  challenge.  We  believe  that  the 
recent  advancements  in  scientific  workstations  will  result  in  a  practical  operational  modeling 
system  much  earlier  than  we  expected. 

C.  SCOPE 

Section  n  discusses  the  tasks  required  to  accomplish  the  project  goals,  and  Section  m 
reviews  the  model  equations,  boundary  conditions,  and  numerical  schemes  used  in  HOTMAC 
and  RAPTAD.  In  Section  IV  we  review  the  MI  diffusion  experiment  and  present  modeling 
results  for  MI87,  MI90,  and  MI91  cases.  In  Section  V,  we  compare  our  model  results  with  those 
obtained  from  current  AF  models  and  RIMPUFF/LINCOM  (Reference  9)  models  developed  at 
RISO  National  Laboratory,  Roskilde,  Denmark. 

Our  experience  with  running  HOTMAC/RAPTAD  on  a  Sun  workstation  is  discussed  in 
Section  VI.  We  studied  the  feasibility  of  incorporating  positive  and  negative  buoyancy  effects 
in  HOTMAC  and  RAPTAD,  and  the  results  are  presented  in  Section  Vn.  Finally,  conclusions 
and  recommendations  for  future  work  are  given  in  Section  VIE. 
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SECTION  n 


TASKS 

A.  PHASE  I 

1.  VAFB  Terrain  Data 

The  contractor  will  enter  digitized  terrain  data/surface  condition  information  for  VAFB 
into  the  HOTMAC/RAPTAD  model.  The  AF  project  officer  will  ensure  that  the  contractor 
receives  the  needed  VAFB  data. 

2.  Model  Runs 

The  contractor  will  operate  HOTMAC/RAPTAD  with  VAFB  terrain  data  for  two  cases 
representing  seasonal  average  meteorology  at  VAFB,  simulating  ground-level  releases  of  a 
neutrally  buoyant  gas.  The  AF  project  officer  will  designate  the  cases  to  be  run. 

3.  Model  Comparisons 

The  contractor  will  compare  output  from  HOTMAC/RAPTAD  for  each  of  the  simu¬ 
lated  ground-level  releases  with  that  of  AFTOX  and  either  OB/DG,  MI,  or  SR.  The  AF  project 
officer  will  designate  which  of  these  models  will  be  used. 

B.  PHASE  n 

1.  HOTMAC/RAPTAD  on  Micro Vax  Computer 

Upon  approval  by  the  AF  project  officer,  the  contractor  will  operate  the  HOT¬ 
MAC/RAPTAD  model  on  a  Micro  Vax  computer  resident  at  LANL.  The  contractor  will  neifonn 
the  model  runs  as  listed  in  paragraphs  in  n.A.2  and  n.A.3. 
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C.  PHASE  m 


The  contractor  will  investigate  incorporating  positive  and  negative  source  buoyancy  effects 
in  RAPTAD  and  determining  the  feasibility  of  incorporating  these  effects. 

Upon  approval  of  the  AF  project  officer,  these  effects  will  be  incorporated  into  the  version 
of  HOTMAC/RAPTAD  resident  at  LANL.  Upon  incorporation  of  these  changes,  the  contractor 
will  operate  HOTMAC/RAPTAD  using  a  test  case  to  determine  the  performance  of  the  model 
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SECTION  m 


MODEL  DESCRIPTION 

Significant  advances  in  engineering  workstation  technologies  have  had  a  phenomenal 
impact  on  computational  fluid  dynamics.  We  have  been  using  a  Sun  4  Workstation  to 
successfully  run  state-of-the-art,  three-dimensional,  atmospheric  transport  and  diffusion  models 
which  were  originally  developed  on  a  supercomputer. 

Our  modeling  system  is  composed  of  two  numerical  codes,  HOTMAC  and  RAPTAD 
(Reference  8).  Only  brief  reviews  of  HOTMAC  and  RAPTAD  are  given  here.  A  detailed 
description  of  the  model  equations,  boundary  conditions,  and  numerical  scheme  is  given  in 
Appendix  A. 

HOTMAC  is  a  mesoscale  atmospheric  model  that  can  forecast  three-dimensional  distri¬ 
butions  of  tvind  speed,  wind  direction,  turbulence,  temperature,  and  water  vapor.  The  basic 
equations  for  HOTMAC  are  the  conservation  equations  for  mass,  momentum,  internal  energy, 
mixing  ratio  of  water  vapor,  and  turbulence  kinetic  energy  (Reference  8). 

Surface  boundary  conditions  are  constructed  from  the  empirical  formulas  by  Dyer  and 
Hicks  (Reference  10)  for  nondimensional  wind  and  temperature  profiles.  The  temperatures  in 
the  soil  layer  are  obtained  by  solution  of  the  heat  conduction  equation.  Appropriate  boundary 
conditions  are  the  heat  energy  balance  at  the  soil  surface  and  specification  of  the  soil  temperature 
at  a  certain  depth.  The  lateral  boundary  values  are  obtained  by  integration  of  the  corresponding 
governing  equations  except  that  variations  in  the  horizontal  directions  are  all  neglected. 

An  initial  wind  profile  at  a  reference  site  in  the  computational  domain  is  first  constructed 
with  the  assumpticHi  of  a  logarithmic  variation  from  the  ground  up  to  the  level  where  the  v/ind 
speed  reaches  an  ambient  value.  Initial  wind  profiles  at  other  grid  locations  are  obtained  from 
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the  winds  determined  above  and  the  conditions  dU/dx  =  0  and  dV/dy  =  0. 

The  governing  equations  are  integrated  by  use  of  the  Alternating  Direction  Implicit  method 
(Reference  11).  A  time  increment  is  chosen  to  be  90  percent  of  the  minimum  value  of  Axi/Ui  , 
where  Ax,  is  a  grid  spacing  and  U;  is  the  velocity  in  the  ith  direction  (Courant-Freidrich-Lewy 
criteria).  To  increase  the  accuracy  of  finite-difference  approximations,  mean  and  turbulence 
variables  are  defined  at  grids  staggered  in  both  the  horizontal  and  vertical  directions.  Mean 
winds,  temperature,  and  water  vapor  vary  most  with  height  near  the  surface.  Nonuniform  grid 
spacings  are  used  in  the  vertical  direction  to  resolve  these  variations  without  introducing  an 
excessive  computational  burden. 

RAPTAD  is  a  Lagrangian  puff  code  based  on  the  Monte  Carlo  statistical  di^sion  process. 
The  center  location  and  standard  deviation  of  concentration  distribution  for  each  puff  are 
computed  by  use  of  wind  and  turbulence  modeled  by  HOTMAC.  Then,  concentration  at  any 
location  is  computed  by  summation  of  concentrations  contributed  by  aU  the  puffs.  RAPTAD  can 
be  used  under  extreme  conditions  with  highly  heterogeneous  wind  and  turbulence  distributions 
where  a  conventional  Gaussian  plume  model  may  fail. 
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SECTION  IV 


MOUNTAIN  IRON  DIFFUSION  DATA  SIMULATIONS 
A.  MOUNTAIN  IRON  DIFFUSION  EXPERIMENTS 

To  test  the  feasibility  of  operating  HOTMACVRAPTAD  with  the  VAFB  terrain  (Tasks 
n.A.l  and  n.A.2),  we  selected  data  from  the  Mountain  Iron  diffusion  experiments  (Reference 
4),  which  were  conducted  at  VAFB  during  1965  and  1966.  The  objective  of  the  MI  program 
was  to  estabhsh  quantitative  diffusion  predictions  for  use  as  range  safety  tools  in  the  “South 
Vandenberg”  (SV)  ballistic  and  space  vehicle  operations  area  (Reference  4).  South  Vandenberg 
is  located  along  the  Cahfomia  coast  approximately  160  km  west- northwest  of  Los  Angeles 
(Figure  1).  The  coastline  is  oriented  in  approximately  a  north-south  direction  along  the  western 
side  of  SV  (Figure  1),  but  changes  abruptly  at  Point  Arguello  to  an  eas;-west  direction.  The 
coastline  gradually  changes  to  north-south  toward  Point  Conception  and  then  changes  again 
to  an  east-west  direction.  The  Santa  Ynez  Mountains  fonn  an  east-west  barrier  along  the 
coastline  far  south  of  SV. 

Another  prominent  terrain  feature  in  the  experimental  area  (Figure  2)  is  Honda  Canyon 
(HC),  walled  by  Target  Ridge  on  the  north  and  Honda  Ridge  on  the  south.  The  distance 
between  the  two  ridges  is  generally  less  than  2.5  km.  The  elevation  difference  between  the 
bottom  of  the  canyon  (less  than  120  meters  above  mean  sea  level,  msl)  and  the  average  height 
of  the  ridges  is  approximately  2(X)  meters. 

Meteorological  masts  10,  50,  1(X),  and  3(X)  feet  high  were  deployed  (Figure  3)  to  measure 
wind  speeds  and  wind  directions.  Temperature  differences  between  6  and  54  feet  above  ground 
level  and  standard  deviations  of  horizontal  wind  directions  were  also  measured  at  the  tracer 
release  points.  The  wind  data  are  used  to  examine  the  performance  of  HOTMAC  in  Section 
rv.B. 
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Figure  1  The  Study  Area,  South  Vandenberg,  is  Located  Approximately  160  km  West 
Northwest  of  Los  Angeles. 


Figtire  2.  Terrain  Features  on  South  Vandenberg. 
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Vertical  profiles  of  wind  speed,  wind  direction,  and  wet  and  dry  temperature  up  to 
2000  meters  above  the  ground  were  measured  by  rawinsondes  at  four  sites,  and  wind  profiles 
were  taken  by  a  wiresonde  at  the  tracer  release  site.  The  soundings  were  taken  hourly  for 
3  to  4  hours  starting  at  1  hour  before  the  tracer  release.  The  vertical  profiles  of  winds  and 
temperatures  are  used  to  initialize  the  model  variables  and  test  the  model  predictions. 

Ruorescent  pigment  zinc  sulfide  particles  with  a  geometric  mean  of  2.5  /zm  were  released 
to  aid  the  understanding  of  transport  and  diffusion  processes  and  the  derivation  of  an  empirical 
formula  for  the  pollutant  concentration  distribution  in  the  SV  area.  The  effective  release  height 
was  2  to  6  meters  above  the  ground.  The  primary  sampler  used  was  a  membrane  filter  inserted 
in  a  disposable  polyethylene  holder.  The  bulk  samples  from  the  field  were  assayed  by  use  of 
the  Rankin  counter,  which  uses  an  alpha  emitter  to  activate  the  fluorescent  pigment  deposited 
on  the  membrane  filter  (Reference  5). 

Accessibility  in  the  rugged  terrain  dictated  the  design  of  the  sampling  grid:  existing  roads 
and  Jeep  trails  provided  a  suitable  and  convenient  network  of  sampling  routes.  The  roads,  in 
general,  follow  creek  banks  and  ridge  tops.  The  chosen  sampler  sites  are  shown  in  Figure 
4.  The  sampling  distances  from  the  primary  source  point  (VIPl)  with  northwest  and  north- 
northwest  prevailing  winds  were  1.5,  2.5,  4,  6,  and  7  km.  The  samplers  along  the  routes  were 
spaced  1/10  mile  (160  meters)  apart  near  the  source,  2/10  mile  (320  meters)  apart  midway,  and 
4/10  mile  (640  meters)  apart  along  the  coastline  far  soutii  of  SV  (Figure  4). 

B.  SIMULATIONS  AND  DISCUSSIONS 

A  total  of  113  tracer  releases  were  made  from  which  8  cases  were  selected  (Reference  12) 
for  detailed  studies  based  on  the  flow-field  categorization  in  the  Plume  Dispersion  Handbook. 

The  wind  data  from  these  eight  cases  were  used  to  evaluate  the  performance  of  the 
diagnostic  wind  model  WOCSS,  >^ds  on  Critical  Streamline  Surfaces  (Reference  12).  More 
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recently.  Rise  National  Laboratory  wind  and  diffusion  models  LINCOM  and  RIMPUFF  were 
evaluated  by  use  of  the  MI87  and  MI90  data  sets  (Reference  13). 

Based  on  these  previous  studies,  we  selected  the  MI87,  MI90,  and  MI91  cases  to  evaluate 
the  performance  of  HOTMAC  and  RAPTAD.  The  following  sections  describe  in  detail  each 
simulation  result  for  the  MI87,  MI90,  and  MI91  cases. 

1.  MI87  Simulation 

The  meteorological  condition  for  the  MI87  case  was  characterized  (Reference  12)  as 
“weak  ambient  and  sea-breeze  flow,  with  little  distinction  between  wind  vectors,  but  sharp 
cooling  with  fog  rolling  in.”  Fog  was  observed  all  morning  at  the  Boat  House  (BH)  site,  but 
none  at  HC  until  after  1500  local  standard  time  (1st).  Fog  progressed  into  the  release  site 
(VIPl),  probably  an  hour  before  the  1310  1st  release  time  (Reference  12). 

The  initial  potential  temperature  profile  was  determined  by  averaging  the  four  upper  air 
soundings  at  1300 1st  at  BH,  Scout  D  (SD),  HC,  and  Building  22  (B22)  (Figure  3).  The  potential 
temperature  lapse  rate  was  approximately  0.044°  C/m  from  the  sea  surface  to  460  meters  above 
msl,  0.0 142°  C/m  between  460  and  960  meters  above  msl,  and  0.0045°  CTm  above  960  metrars 
above  msl.  )Mnd  speed  and  direction  were  determined  by  examining  five  upper  air  soundings 
(four  locations  mentioned  above  plus  VIPl)  at  1300 1st  Initial  upper  air  wind  speed  and  wind 
direction  above  the  boundary  layer  were  estimated  to  be  2  m/s  and  225  degrees,  respectively. 

The  computational  domain  is  40  x  48  km^  with  a  horizontal  grid  spacing  of  1  km. 
To  resolve  the  details  of  topography  in  the  vicinity  of  the  release  site,  we  decided  to  nest  a 
fine-resolution  grid  15  x  16  km^  with  a  horizontal  grid  spacing  of  0.5  km. 

Integration  started  at  0500  1st,  June  13,  1966,  and  continued  for  over  12  hours.  The 
plume  was  released  at  1310  1st  for  30  minutes  as  was  done  in  the  experiment  The  plume 
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was  followed  for  4  hours  in  the  model  computation.  By  that  time,  the  plume  had  moved  far 
from  the  sampling  areas. 

Figure  5  shows  the  modeled  horizontal  wind  vectors  in  the  inner  computational  grid 
at  10  meters  above  the  ground  at  1300  1st,  June  13  (Julian  date  164).  Although  the  upper 
air  wind  direction  is  225  degrees  (southwest),  upslope  flows  develop  in  the  surface  layer  as  a 
result  of  heating  at  the  sloped  surfaces. 

The  modeled  wind  distribution  (Figure  5)  agrees  with  the  observation  (Figure  6).  The 
observed  winds  show  much  more  variation  in  space  than  the  modeled  winds.  Observations 
near  each  other  show  considerable  variation  in  direction  and  magnimde.  On  the  other  hand,  the 
modeled  wind  field  varies  more  slowly  in  space  than  the  observed  wind  field  since  the  model 
neglects  subgrid-scale  variations  of  the  surface  (the  grid  resolution  is  500  meters).  Nevertheless, 
the  simulations  successfully  reproduced  many  observed  features. 

Less  satisfactory  results  are  obtained  m  comparison  of  the  vertical  profiles  of  the 
modeled  wind  speed  and  wind  direction  with  observations  (Figure  7).  Wind  speed  and  wind 
direction  vary  greatly  in  space  and  time  when  the  wind  speed  is  small.  The  observations  are 
instantaneous  values  whereas  the  modeled  results  are  ensemble  averages.  On  the  other  hand, 
potential  temperature  profiles  are  relatively  stationary  unless  synoptic-scale  disturbances  such 
as  fronts  pass  through  the  measurement  area. 

Significant  changes  in  the  modeled  wind  direction  occurred  at  around  600  meters  above 
the  ground.  These  changes  are  caused  by  the  mass  conservation  constraint  to  compensate  for 
the  divergence  and  convergence  of  the  wind  distributions  in  the  boundary  layer  (Figure  5). 
Observed  wind  direction  profiles  appear  to  support  such  variations,  but  the  changes  appear  to 
occur  at  heights  much  closer  to  the  ground  than  those  modeled. 
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Figure  5.  Modeled  Horizontal  Wind  Vectors  for  the  MI87  Cetse  at  10  meters  above  the 
Ground  at  1300 1st,  June  13, 1966.  Wind  vectors  at  every  other  grid  point  are 
plotted.  Terrain  is  contoured  by  solid  lines  with  an  increment  of  200  meters. 
Dashed  lines  indicate  contours  halfway  between  the  solid  contours. 
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Figure  6.  Same  as  in  Figure  5  except  Observed  Wind  Vectors  in  the  Surface  Layer  are  Shown. 
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Figure  7.  Vertical  Profiles  of  the  Modeled  Wind  Speed,  Direction,  and  Potential  Tem¬ 
perature  for  the  MI87  Case  at  1300  1st,  June  13,  1966  at  (a)  BH,  (b)  HC, 
(c)  SD,  (d)  VIPl,  and  (e)  B22.  Solid  circles  indicate  observations. 
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Figure  7.  Vertical  Profiles  of  the  Modeled  Wind  Speed,  Direction,  and  Potential  Tem¬ 
perature  for  the  MI87  Case  at  1300  1st,  Jime  13,  1966  at  (a)  BH,  (b)  HC, 
(c)  SD,  (d)  VlPl,  and  (e)  B22.  Solid  circles  indicate  observations. 
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Figure  7.  Vertical  Profiles  of  the  Modeled  Wind  Speed,  Direction,  and  Potential  Tem¬ 
perature  for  the  MI87  Case  at  1300  1st,  June  13,  1966  at  (a)  BH,  (b)  HC, 
(c)  SD,  (d)  VIPl,  and  (e)  B22.  Solid  circles  indicate  observations. 
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Figure  7.  Vertical  Profiles  of  the  Modeled  Wiod  Si>eed,  Direction,  and  Potential  Tem¬ 
perature  for  the  MI87  Case  at  1300  kt,  June  13,  1966  at  (a)  BH,  (b)  HC, 
(c)  SD,  (d)  VTPl,  and  (e)  B22.  Solid  circles  indicate  observations. 
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Figure  7.  Vertical  Profiles  of  the  Modeled  Wind  Speed,  Direction,  and  Potential  Tem- 
peratiure  for  the  MI87  Case  at  1300  1st,  June  13,  1966  at  (a)  BH,  (b)  HC, 
(c)  SD,  (d)  VIPl,  and  (e)  B22.  Solid  circles  indicate  observations. 
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Figure  8  shows  the  modeled  ground-level  exposure  contours  and  Figure  9  shows  the 
corresponding  observation.  An  exposure  value  is  computed  by  integrating  concentration  with 
time  until  the  change  in  concentration  becomes  negligible.  Then  it  is  divided  by  the  total 
amount  of  the  materials  released.  The  unit  of  exposure  value  is  s/m^.  Figure  9  also  shows  the 
observed  wind  speeds  and  wind  directions  at  the  ground  stations.  The  observed  wind  direction 
close  to  the  release  site  is  west-southwesterly,  but  it  changes  to  westerly  at  the  station  slightly 
north  of  the  release  site.  Our  simulation  (Figure  5)  indicates  that  the  wind  direction  is  close 
to  westerly  at  the  release  site.  The  observed  plume  apparently  transported  to  the  east-northeast 
direction  despite  the  fact  ±at  wind  directions  measured  at  the  ground  stations  suggest  the  plume 
should  be  transported  to  the  east-southeast,  which  is  the  case  for  the  modeled  plume. 

Although  the  modeled  plume  direction  did  not  match  the  observed  one,  the  modeled 
ground-level  concentrations  along  the  plume  axis  are  in  good  agreement  with  observations,  as 
shown  in  Table  1.  It  is  not  known  why  the  observation  at  3310  meters  from  the  source  shows 
the  largest  value  among  the  observations. 

TABLE  1:  COMPARISON  BETWEEN  THE  MODELED  AND 
OBSERVED  EXPOSURES  FOR  MI87  CASE 


Modeled 

Observed 

Distance  from 
the  source 
(m) 

Exposure 

(s/m^) 

Distance  from 
the  source 
(m) 

Exposure 

(s/m^) 

716 

8.04  X  10-^ 

720 

3.97  X  10-^ 

1253 

2.55  X  10-^ 

1260 

1.08  X  10-^ 

3312 

9.6  X  10-'^ 

3310 

5.39  X  10-<^ 

4403 

7.8  X  10  '^ 

4400 

3.02  X  lO-"^ 
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Figure  8.  Modeled  Ground  Level  Exposure 
Case.  Units  are  s/m®. 


Finally,  the  modeled  horizontal  wind  vectors  in  the  outer  grid  at  14  meters  above  the 
ground  at  1300  1st  are  shown  in  Figure  10  which  indicates  that  a  sea-breeze  front  penetrates 
deep  inland  since  the  sea  breeze  is  enhanced  by  the  flows  developed  over  the  slopes  located 
near  the  eastern  boundary  of  the  study  area. 

2.  MI90  Simulation 

Meteorological  conditions  for  the  MI90  case  are  significantly  different  from  those  for 
the  MI87  case.  First,  the  prevailing  wind  speed  for  MI90  was  close  to  13.5  m/s  compared  with 
2  m/s  for  MI87.  Second,  for  MI90,  a  deep,  well-mixed  layer  existed  for  the  first  500  meters 
above  the  ground,  which  probably  resulted  from  strong  mixing  caused  by  turbulence  generated 
by  high  winds.  On  the  other  hand,  the  potential  temperature  profiles  for  MI87  exhibited  shallow 
(100  meters )  mixed-layer  depths,  typical  on  the  west  coast  as  a  result  of  large-scale  subsidence. 

The  initial  potential  temperature  profile  was  determined  firom  the  sounding  at  B22  site. 
The  potential  temperature  lapse  rate  was  almost  neutral  (0.0(X)2°C/m)  from  the  surface  to  670 
meters  above  msl,  0.061  l°Cym  between  670  and  880  meters  above  msl,  and  0.(X)70°C/m  above 
880  meters  above  msl.  Initial  wind  speed  and  direction  were  determined  by  examination  of 
three  upper  air  soundings  (HC,  SD,  and  B22)  at  2300  1st.  The  wind  profile  at  BH  showed  an 
extremely  high  wind  speed  (25  m/s)  compared  with  those  (13  m/s)  at  other  sites.  Thus,  the 
BH  profile  was  not  considered  in  the  determination  of  initial  wind  speed  of  13.5  m/s  and  wind 
direction  of  332  degrees  which  were  used  in  the  model 

The  initial  relative  humidity  profile  was  determined  from  the  composite  profiles  at 
four  sites  (BH,  HC,  SD,  and  B22).  The  relative  humidity  was  80  percent  at  the  ground  and 
increased  linearly  with  height  to  92  percent  at  670  meters  above  msl.  Then,  it  decreased  rapidly 
to  38  percent  at  880  meters  above  msl.  A  constant  value  (38  percent)  was  used  in  the  layers 
above  880  meters  above  msl 
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Figure  10.  Same  as  in  Figure  5  except  in  the  Larger  Computational  Domain. 
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The  computational  domain  is  exacdy  the  same  as  for  the  MI87  case:  40  x  48  km^ 
with  a  horizontal  grid  spacing  of  1  km,  and  the  corresponding  values  for  the  nested  inner  grid 
are  15  X  16  km  ^  and  0.5  km,  respectively. 

Integration  started  at  0600  1st,  June  21,  1966,  and  continued  for  24  hours  until  0600 
1st  the  following  day.  The  plume  was  released  for  30  minutes  starting  at  2300  1st,  June  21, 
1966.  The  plume  was  sampled  for  2  hours  in  the  model  computation.  By  the  end  of  2  hours  of 
sampling,  the  plume  had  left  the  computational  domain,  and  the  surface  concentration  values 
were  practically  zero. 

The  modeled  horizontal  wind  vectors  (Figure  11)  at  10  meters  above  the  ground  at 
2300,  June  21,  1966  (Julian  date  172),  show  very  little  variation  because  the  prevailing  wind 
is  strong  enough  (13.5  m/s)  to  eliminate  differential  forcing  generated  by  the  topography.  The 
observed  horizontal  wind  vectors  (Figure  12),  on  the  other  hand,  show  considerable  spatial 
variations  partly  because  the  measurement  heights  at  each  site  are  different 

Vertical  profiles  of  the  modeled  wind  speed,  wind  direction,  and  potential  temperature 
(Figures  13)  are  in  fair  agreement  with  observations.  The  observed  wind  speed  profiles  at  HC 
and  BH  are  quite  different  fiom  those  at  B22  and  VlPl,  which  exhibit  normal  variations.  We 
are  unable  to  find  physical  reasons  for  the  extraordinary  wind  profiles  at  HC  and  BH,  which 
were  consistent  over  the  time  period  between  2200  1st,  June  21,  and  0300  1st,  June  22,  1966. 

Wfith  the  wind  and  turbulence  distributions  computed  by  HOTMAC,  RAPTAD  simu¬ 
lated  the  ground-level  concentration  contours  (Hgure  14),  which  are  in  good  agreement  with 
the  observations  (Figure  15).  The  observed  exposures  at  1870  and  5620  meters  from  the  source 
were  5.08  x  10~®  and  4.97  x  10“^  s/m^,  respectively.  The  cotiesponding  modeled  values  at 
1895  and  5603  meters  from  the  source  are  1.83  x  10“®  and  3.45  x  10”^  s/m^,  respectively. 
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Figure  11.  Modeled  Horizontal  Wind  Vectors  for  the  MI90  Case  at  10  meters  above  the 
Groimd  at  2300 1st,  Jime  21, 1966.  Wind  vectors  at  every  other  grid  points  are 
plotted.  Terr^  is  contoured  by  solid  lines  with  an  increment  of  200  meters. 
Daahed  lines  indicated  contours  halfway  between  the  solid  contours. 
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Figure  13.  Vertical  Profiles  of  the  Modeled  Wind  Speed,  Direction,  and  Potential  Tem- 
peratiue  for  the  MI90  Case  at  2300  1st,  June  21,  1966  at  (a)  BH,  (b)  HC, 
(c)  SD,  (d)  VIPl  and  (c)  B22.  Solid  circles  indicate  observations. 
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Figure  13.  Vertical  Profiles  of  the  Modeled  Wind  Speed,  Direction,  and  Potential  Tem¬ 
perature  for  the  MI90  Case  at  2300  1st,  June  21,  1966  at  (a)  BH,  (b)  HC, 
(c)  SD,  (d)  VIPl  and  (e)  B22.  Solid  circles  indicate  observations. 
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Figiire  13.  Vertical  Profiles  of  the  Modeled  Wind  Speed,  Direction,  and  Potential  Tem¬ 
perature  for  the  MI90  Case  at  2300  1st,  June  21,  1966  at  (a)  BH,  (b)  HC, 
(c)  SD,  (d)  VIPl  and  (e)  B22.  Solid  circles  indicate  observations. 
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Figure  13.  Vertical  Profiles  of  the  Modeled  Wind  Speed,  Direction,  and  Potential  Tem¬ 
perature  for  the  MI90  Case  at  2300  bt,  June  21,  1966  at  (a)  BH,  (b)  HC, 
(c)  SD,  (d)  VIPl  and  (e)  B22.  Solid  circles  indicate  observations. 
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Figtire  13  Vertical  Profiles  of  the  Modeled  Wind  Speed,  Direction,  and  Potential  Tem¬ 
perature  for  the  MI90  Case  at  2300  1st,  June  21,  1966  at  (a)  BH,  (b)  HC, 
(c)  SD,  (d)  VIPl  and  (e)  B22.  Solid  circles  indicate  observations. 
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3.  MI91  Simulation 


The  MI91  release  was  at  0203 1st,  June  22,  1966  which  is  only  3  hours  after  the  MI90 
release.  However,  the  wind  direction  changed  considerably  to  northerly,  as  much  as  30  degrees 
from  the  MI90  case.  Based  on  the  MI91  data  log,  the  wind  direction  veering  from  northwest 
to  north  occurred  around  midnight  June  21,  1966  and  was  more  evident  in  the  layers  between 
300  and  1300  meters  above  msl.  At  the  present  time,  there  is  no  valid  documentation  to 
explain  how  the  sudden  wind  direction  shift  occurred  except  that  the  MI90  data  log  suggested 
a  possible  postfrontal  (i.e.,  synoptically  forced)  effect. 

The  maximum  wind  veering  occurred  at  around  800  meters  above  msl,  and  a  sine 
curve  appears  to  represent  well  the  vertical  profiles  of  the  wind  direction  change  between  300 
and  1300  meters  above  msl.  Thus,  wind  directions  at  the  grid  levels  between  300  and  1300 
meters  above  msl  were  nudged  toward  the  observed  wind  directions.  We  started  the  HOTMAC 
computation  with  the  wind  direction  nudging.  The  output  at  2300  from  the  MI90  simulation 
was  used  to  initialize  the  model  variables. 

The  modeled  horizontal  wind  vectOTS  (Figure  16)  clearly  indicate  that  the  wind 
directions  have  shifted  from  the  corresponding  values  for  the  MI90  case  (Figure  11).  The 
nudging  was  not  applied  to  the  layers  in  the  first  300  meters  above  the  ground,  but  the  surface 
layer  winds  (Figure  16)  are  clearly  influenced  by  the  nudging.  This  good  communication 
between  the  surface  and  upper  level  winds  apparently  resulted  from  the  strong  vertical  mixing 
resulting  from  extraordinary  large  turbulence  (Figure  17)  for  a  nocturnal  period.  The  high 
prevailing  wind  speed  (13.5  m/s)  is  responsible  for  the  large  turbulence  values. 

Unlike  the  modeled  horizontal  wind  vectors  at  10  meters  above  the  ground  (Figure 
16),  the  observed  surface  wind  distribution  for  MI91  (Figure  18)  is  not  much  different  from 
the  corresponding  distribution  for  MI90  (Figure  12). 
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Figvire  1C.  Modeled  Horizontal  Wind  Vectors  for  the  MI91  Case  at  10  meters  above  the 
Ground  at  0200  1st,  Jxine  22,  1966.  Wind  vectors  are  plotted  at  every  other 
grid  point.  Terrain  is  contoured  by  sohd  hnes  with  an  increment  of  200  meters. 
Dashed  lines  indicate  contours  halfway  between  the  soUd  contours. 
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Figure  17.  Vertical  Profiles  of  the  Modeled  Standard  Deviation  of  Vertical  Wind  Compo¬ 
nent,  Turbulence  Kinetic  Energy,  and  Standard  Deviation  of  Potential  Tem- 
peratine  at  VIPl  at  0200  1st,  June  22,  1966. 
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Figure  18.  Same  as  in  Figure  16  except  Observed  Wind  Vectors  in  the  Surface  Layer  are 
Shown. 
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Vertical  profiles  of  the  modeled  wind  speed,  wind  direction,  and  potential  temperatiue 
(Figure  19)  are  in  fair  agreement  with  observations:  the  modeled  wind  direction  with  the 
nudging  are  in  good  agreement  with  measurements,  but  wind  speed  and  temperature  profiles  are 
quite  different  from  the  observations.  The  MI91  data  log  hinted  that  the  postfrontal  disturbance 
might  have  caused  the  veering  and  strong  warming  in  the  boundary  layer,  but  the  description 
was  not  good  enough  for  us  to  implement  the  large-scale  variations  in  HOTMAC. 

Using  the  wind  and  turbulence  distributions  computed  by  HOTMAC,  RAPTAD  sim¬ 
ulated  ground-level  concentration  contours  (Figure  20)  which  are  in  good  agreement  with  the 
observations  (Figure  21).  The  observations  appear  to  indicate  a  minimum  value  in  Honda 
Canyon,  and  the  plume  axis  is  more  northerly  than  the  modeled  plume  axis.  The  modeled  ex¬ 
posures  along  the  plume  axis  are  in  good  agreement  with  the  observations  at  nearby  locations 
as  shown  in  Table  2. 


TABLE  2:  COMPARISON  BETWEEN  THE  MODELED  AND  OBSERVED 
EXPOSURES  FOR  MI91  CASE. 

Modeled  Observed 


Distance  from 
the  soiuce 
(m) 

Exposure 

(s/m^) 

Distance  from 
the  source 
(m) 

Exposure 

(s/m^) 

1254 

3.30  X  10-^ 

1270 

4.09  X  10-^ 

2672 

1.03  X  10-^ 

2670 

1.54  X  10-® 

4253 

6.98  X  lO  "^ 

4220 

9.11  X  10-^ 

5005 

5.45  X  10-'' 

4950 

2.42  X  10-^ 

7139 

3.15  X  lO-’^ 

7130 

2.02  X  lO  ”^ 

10030 

1.06  X  lO  "^ 

10000 

5.07  X  10-* 
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Figure  19.  Vertical  Profiles  of  the  Modeled  Wind  Speed,  Direction,  and  Potential  Tem¬ 
perature  for  the  MI91  Case  at  0200  1st,  June  17,  1966  at  (a)  BH,  (b)  HC, 
(c)  SD,  (d)  VlPl,  and  (e)  B22.  Solid  circles  indicate  observations. 
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Figure  19.  Vertical  Profiles  of  the  Modeled  Wind  Speed,  Direction,  and  Potential  Tem¬ 
perature  for  the  MI91  Case  at  0200  1st,  June  17,  1966  at  (a)  BH,  (b)  HC, 
(c)  SD,  (d)  VIPl,  and  (e)  B22.  Solid  circles  indicate  observations. 
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Figure  19.  Vertical  Profiles  of  the  Modeled  Wind  Speed,  Direction,  and  Potential  Tem¬ 
perature  for  the  MI91  Case  at  0200  1st,  June  17,  1966  at  (a)  BH,  (b)  HC, 
(c)  SD,  (d)  VIPl,  and  (e)  B22.  Solid  circles  indicate  observations. 
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Figiire  19.  Vertical  Profiles  of  the  Modeled  Wind  Speed,  Direction,  and  Potential  Tem¬ 
perature  for  the  MI91  Case  at  0200  1st,  June  17,  1966  at  (a)  BH,  (b)  HC, 
(c)  SD,  (d)  VIPl,  and  (e)‘B22.  Solid  circles  indicate  observations. 
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Figure  19.  Vertical  Profiles  of  the  Modeled  Wind  Speed.  Direction,  and  Potential  Tem- 
p>erature  for  the  MI91  Case  at  0200  1st,  June  17,  1966  at  (a)  BH,  (b)  HC, 
(c)  SD,  (d)  VIPl,  and  (e)  B22.  Sohd  circles  indicate  observations. 
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The  HOTMAC  and  RAPTAD  modeling  system  produced  surface  exposures  along  the 
plume  axis  at  least  as  good  as  those  obtained  by  diagnostic  models  and  an  empirical  formula 
specifically  adjusted  to  the  MI  data.  The  real  strength  of  a  prognostic  model  such  as  HOTMAC 
is  realized  where  wind  data  are  not  available.  Under  such  conditions,  only  prognostic  models 
can  provide  a  realistic  solution  which  considers  complex  surface  conditions.  In  the  following 
section,  the  HOTMAC  and  RAPTAD  results  are  compared  with  those  obtained  by  diagnostic 
models  and  AF  models  currently  used  at  VAFB. 
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SECTION  V 


MODEL  COMPARISONS 

There  were  a  few  attempts  in  the  past  to  simulate  the  MI  wind  and  concentration  data. 
Hunter  (Reference  12)  investigated  the  performance  of  a  diagnostic  wind  model  WOCSS  for 
eight  cases  of  the  MI  data.  Thykier — Nielsen  et  al.  (Reference  13)  tested  a  diagnostic  wind 
model  (LINCOM)  and  a  puff  model  (RIMPUFF)  for  two  cases  of  the  MI  data.  The  data 
common  to  both  studies  were  the  MI87  and  MI90  cases.  Therefore,  we  use  these  two  cases  to 
compare  the  performances  of  WOCSS,  LINCOM,  RIMPUFF,  HOTMAC,  and  RAPTAD,  and 
discuss  the  strengths  and  weaknesses  of  each  model.  The  present  comparisons  are  qualitative 
and  informal  because  none  of  the  model  results  presented  here  has  experienced  the  formal 
review  normally  required  for  journal  publication. 

A.  WIND  MODEL 

There  are  two  classes  of  wind  models:  diagnostic  and  prognostic.  Diagnostic  vind  models 
are  based  on  measurements  and  interpolation.  Most  of  these  models  adjust  initially  interpolated 
wind  fields  to  satisfy  a  mass  conservation  equation.  The  accuracy  of  the  model  results  is  greatly 
dependent  upon  the  quality  and  representativeness  of  measurements. 

Diagnostic  wind  models  have  been  used  almost  exclusively  for  emergency  response 
management  because  they  are  simple  and  fast.  The  critical  weakness  of  this  class  of  models, 
however,  is  the  lack  of  predictive  capability.  Plume  models  based  on  a  diagnostic  wind  model 
must  assume  that  the  current  wind  field  will  persist  until  it  is  updated  by  future  measurements. 
The  assumption  of  persistency  is  adequate  if  the  measurements  are  made  frequently  and  if  no 
abrupt  change  in  wind  field  is  expected.  However,  such  conditions  limit  the  area  of  concern 
to  the  immediate  vicinity  of  the  source  and  exclude  predictions  during  the  transitional  periods. 
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The  question  of  the  representativeness  of  measurements  further  limits  the  applications 
of  diagnostic  wind  models  to  simple  terrains,  because  complex  surface  boundaries  require  a 
prohibitively  large  number  of  wind  measurements. 

Prognostic  models,  on  the  other  hand,  solve  a  set  of  time-dependent  physical  equations, 
such  as  conservation  equations  of  momentum,  internal  energy,  and  mixing  ratio  of  water  vapor. 
Prognostic  models  forecast  three-dimensional  distributions  of  wind  field  that  become  the  input 
to  transport  and  diffusion  models. 

At  present,  prognostic  models  are  not  used  in  emergency  response  systems  because  a  large 
amount  of  computation  time  is  required.  The  situation  should  be  improved  considerably  in  the 
near  future,  because  the  capabilities  of  engineering  workstations  are  expected  to  advance  at 
an  astonishing  rate. 

Figure  22  shows  the  modeled  horizontal  wind  vectors  for  the  MI87  case  at  4  meters  above 
ground  level  (agl)  at  1310  1st  by  using  WOCSS  diagnostic  wind  model  (Reference  12).  An 
array  of  50  x  80  x  6  with  500  meter  grid  spacing  and  the  terrain  data  of  a  500  meter  resolution 
were  used.  )^d  vectors  at  every  fourth  grid  point  are  shown  in  Figure  22. 

Figure  23  shows  the  modeled  wind  vectors  for  the  MI87  case  at  10  meters  agl  at  1310  1st 
by  using  LINCXDM  (Reference  13),  which  is  a  linear,  diagnostic,  spectral,  and  potential-flow 
wind  model. 

Both  WOCSS  and  LINCOM  used  the  same  tower  data  and  wind  fields  were  adjusted  to 
satisfy  mass  consistency.  Indeed,  wind  distributions  shown  in  Figures  22  and  23  are  in  close 
agreement  with  each  other,  although  LINCOM  (Hgure  23)  produced  slightly  more  spatial 
variations  than  WOCSS  (Figure  22). 

As  mentioned  earlier,  a  diagnostic  wind  model  is  no  more  than  a  scheme  for  interpolation 
of  observed  winds  at  specified  grid  points.  For  this  reason,  the  spatial  variations  of  wind 
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Figure  22.  Modeled  Horizontal  Wind  Vectors  for  the  MI87  Case  at  4  meters  above  Groimd 
Level  by  Using  WOCSS  Diagnostic  Wind  Model.  Meteorological  measurement 
stations  are  identified  by  numbers. 


Figuie  23.  Modeled  Horizontal  Wind  Vectors  for  the  MI87  Cast  at  10  meters  above 
Ground  Level  by  Using  LINCOM  Diagnostic  Wind  Model. 
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vectors  in  Figures  22  and  23  are  confined  mainly  to  the  regions  where  data  are  available.  Mass 
consistency  and  other  physical  constraints  used  in  diagnostic  wind  models  can  affect  the  wind 
vectors  away  from  the  measurement  sites. 

A  prognostic  wind  model,  on  the  other  hand,  is  based  on  a  set  of  physical  equations 
and  does  not  require  the  measurement  of  surface  winds  for  forecast  Observations  are  used 
to  verify  the  predictions. 

The  prognostic  model  results  (Figure  5)  generally  agree  with  the  diagnostic  model  results 
(Figures  22  and  23)  in  the  area  where  wind  measurements  are  available,  ^^^d  vectors  in  the 
prognostic  model  are  determined  from  the  balance  among  the  pressure  force,  turbulent  mixing, 
and  Coriolis  force.  Therefore,  prognostic  models  can  provide  wind  distributions  in  regions 
where  no  measurements  are  available.  Our  wind  vectors  (Figure  10)  in  the  outer  grid  are 
quite  different  from  the  wind  fields  computed  by  WOCSS  (Figure  22)  and  LINCOM  (Figure 
23)  in  the  areas  where  topographic  features  are  prominent.  The  diagnostic  models  do  not 
incorporate  the  pressure  forces  generated  by  the  sloped  surfaces.  Thraefore,  the  variations  of 
wind  distributions  based  on  a  diagnostic  wind  model  must  come  from  the  wind  data  that  reflect 
the  surface  inhomogeneity. 

Figures  24  and  25  are  the  wind  vectors  for  the  MI90  case  modeled  by  WOCSS  and 
LINCOM,  respectively.  The  wind  directions  are  uniform  (northwesterly)  except  in  the  areas 
where  tower  wind  data  were  available. 

Wind  vectors  simulated  by  HOTMAC  (Figure  11)  show  very  little  variations  in  space 
because  the  strong  (13.5  m/s)  prevailing  wind  overpowered  the  perturbations  generated  by 
differential  cooling  over  the  sloped  surface.  Figure  11  is  quite  different  from  Figure  5,  in 
which  the  prevailing  wind  was  small  (2  m/s)  and  the  modeled  surface  wind  vectors  were 
perturbed  significantly  by  differential  heating  over  the  sloped  surfaces.  The  comparison  between 
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Figure  24.  Same  as  in  Figure  22  except  for  the  MI90  Case. 
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Figure  25.  Same  as  in  Figure  23  except  for  the  MI90  Case. 
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Figures  5  and  11  demonstrate  the  HOTMAC’s  capabilities  in  realistically  incorporating  the 
effects  of  large-scale  pressure  force  and  local-scale  pressure  force  due  to  differential  heating. 
B.  DIFFUSION  MODEL 

Thykier — Nielsen  et  al.  (Reference  13)  applied  a  diffusion  model  (RIMPUFF)  to  simulate 
the  fluorescent  particle  concentration  distributions  at  the  ground  for  the  MI87  and  MI90  cases. 
RIMPUFF  is  a  Lagrangian  puff  model  that  simulates  time-variable  continuous  releases  by 
sequentially  releasing  a  series  of  Gaussian-shaped  puffs.  Mean  wind  fields  were  provided  by 
LINCOM  and  turbulence  velocities  were  based  on  the  local  energy  of  the  fluctuating  horizontal 
wind  components.  The  turbulent  energies  were  estimated  from  the  climatological  study  of  flow 
and  turbulence  at  VAFB  carried  out  in  cormection  with  handbook-study  (Reference  13).  The 
turbulence  velocity  was  added  to  the  mean  value  to  determine  a  puff  velocity  at  the  center  of 
mass.  The  procedure  is  similar  to  that  used  in  RAPTAD  except  that  turbulence  velocities  in 
RAPTAD  are  provided  by  HOTMAC. 

Figures  26  and  27  show  the  concentration  distributions  of  fluorescent  particles  at  ground 
level  as  modeled  by  RIMPUFF  for  the  MI87  and  MI90  cases,  respectively.  Thykier — ^Nielsen  et 
al.  (Reference  13)  do  not  explain  the  values  for  the  concentration  contours,  nor  do  they  present 
a  comparison  of  the  modeled  and  observed  concentration  values.  However,  Thykier — ^Nielsen 
et  al.  did  comment  on  their  results,  saying,  “For  all  of  the  eight  tracer  experiments  studied,  our 
simulated  mean-plume  spread  resembles  well  in  form  and  shape  the  corresponding  measured 
tracer  concentrations,  at  least  for  the  relatively  short  range  where  concentration  measurements 
were  detectable  during  the  Motmtain  Itot  tracer  experiments.” 

Finally,  we  compare  the  performance  of  HOTMAC  and  RAPTAD  with  that  of  the 
MI  model  currently  used  at  VAFB.  The  MI  model  is  an  empirical  equation  that  relates  a 
surface  concentration  value  to  distance,  standard  deviation  of  horizontal  wind  fluctuation,  mean 
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Figure  27.  Same  as  in  Figure  26  except  for  the  MI90  Case. 
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velocity,  and  temperature  gradient  in  the  surface  layer.  Empirical  coefficients  were  determined 
to  fit  best  the  MI  diffusion  data. 

A  simple  distance-dependent  formula  resulted  in  good  fit  to  the  observations,  with  63 
percent  within  a  factor  of  2  of  the  mean  and  92  percent  within  a  factor  of  4.  We  applied  the 
same  formula  at  the  sampling  sites  for  MI87,  MI90,  and  MI91  cases  and  obtained  exposure 
values.  The  results  are  compared  with  observations  in  Figure  28,  where  HOTMAC/RAPTAD 
predictions  are  also  plotted.  The  predicted  values  are  within  a  factor  of  2  of  the  observations 
between  the  inner  two  dashed  lines  and  within  a  factor  of  4  between  the  outer  two  dashed 
lines.  Two  data  points  outside  of  the  factor  of  4  boundary  resulted  from  the  unreasonably  large 
observed  exposure  value  reported  for  the  MI87  case. 

Figure  28  indicates  that  HOTMAC  and  RAPTAD  predictions  are  at  least  as  good  as  those 
made  by  the  MI  model  which  was  specifically  tuned  with  the  MI  diffusion  data. 
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Figiire  28.  Scatter  Diagram  of  the  Predicted  and  Observed  Fluorescent  Particle  Concen¬ 
tration  (nondimensional)  at  Ground  Level.  The  inner  dashed  lines  indicate  a 
factor  of  2  and  the  outer  dashed  hnes  indciate  a  factor  of  4. 


SECTION  VI 


APPLICATIONS  ON  DESKTOP  COMPUTERS 

Recent  advances  in  desktop  computer  capabilities,  particularly  those  of  an  engineering 
workstation,  are  astonishing.  A  high-performance  workstation  has  reportedly  exceeded  a  super¬ 
computer  in  certain  scalar  operations.  The  affordabUity  and  portability  of  the  desktop  computer 
have  opened  the  door  to  many  applications  that  were  previously  considered  impossible. 

One  logical  area  for  such  applications  is  in  upgrading  the  toxic  hazard  modeling  capabilities 
for  emergency  response  management  Scientists  at  Los  Alamos  National  Laboratory  are  in  the 
process  of  developing  three-dimensional  forecasting  models  that  run  on  a  desktop  computer. 

In  this  section  we  will  discuss  our  experience  in  running  HOTMAC  and  RAPTAD  on  a 
Sun  Microsystems  workstation  and  a  MicroVax  computer. 

A.  SUN  4/110  WORKSTATION 

The  HOTMAC  code  for  a  desktop  computer  is  slightly  different  from  that  for  a  super¬ 
computer.  First,  a  single  grid  is  used.  The  number  of  cells  is  generally  smaller.  This  requires 
either  a  coarser  terrain  resolution  or  a  smaller  domain,  depending  on  the  simulation.  Finally, 
the  vertical  profile  of  the  initial  potential  temperature  may  be  represented  by  two  segments, 
whereas  three  are  used  in  the  supercomputer  version. 

We  ran  HOTMAC  and  RAPTAD  for  the  MI87  case  on  a  Sun  4/110  workstation.  The 
computational  domain  matches  the  outer  grid  used  in  Section  IV,  i.e.,  40  x  48  km^.  The 
horizontal  grid  spacing  was  2  km  as  opposed  to  I  km  in  the  Cray  version.  Potential  temperature 
lapse  rates  of  0.04°  C/m  from  the  sea  surface  to  460  meters  above  msl  and  0.0044°  C/m  above 
460  meters  above  msl  were  used.  In  the  Cray  version  an  additional  lapse  rate  of  0.0142°C/m 
was  used  from  460  to  960  meters  above  msl.  Other  meteorological  conditions  were  identical. 
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Figure  29  shows  the  modeled  horizontal  wind  vectors  14  meters  above  the  ground  at 
0700, 1400,  and  2300 1st  June  14, 1966.  Although  the  upper  air  wind  direction  was  225  degrees 
(southwesterly),  downslope  flows  and  land  breeze  developed  in  the  surface  layer  due  to  cooling 
at  the  ground.  Sumise  on  June  14,  was  shortly  before  0500  1st. 

As  the  sun  warms  the  ground,  the  air  temperature  in  the  surface  layer  quickly  increases  as 
a  result  of  turbulent  mixing.  The  air  temperature  over  a  sloping  surface  becomes  higher  than 
the  temperature  at  the  same  elevation  but  away  from  the  ground.  The  temperature  difference 
results  in  a  horizontal  pressure  gradient:  3.  low  pressure  is  normally  located  at  high  elevation. 
This  pressure  force  produces  upslope  flows  and  sea  breezes  that  converge  over  ridges  as  seen 
in  the  horizontal  wind  vector  distribution  at  1400  1st 

The  ground  starts  cooling  shortly  before  sunset  (around  2000  1st)  and  the  air  temperature 
close  to  the  ground  over  a  sloped  surface  falls  below  the  air  temperature  at  the  same  elevation 
but  away  from  the  surface.  Thus,  the  direction  of  the  pressure  gradient  at  night  is  opposite 
to  its  daytime  counterpart.  However,  the  magnitude  of  the  pressure  gradient  force  at  night  is 
normally  much  smaller  than  in  the  daytime,  particularly  in  summer,  for  the  following  reasons. 

First,  turbulence  intensity  in  the  surface  layer  during  the  nocturnal  period  is  extremely 
small  compared  with  the  daytime  value.  Therefore,  only  the  air  relatively  close  to  the  ground 
is  cooled  by  turbulent  mixing.  The  depth  of  the  nocturnal  stable  boundary  layer  is  on  the  order 
of  100  meters,  much  smaller  than  that  of  the  daytime  convective  boundary  layer,  which  is  at 
least  several  hundred  meters  high. 

Second,  the  solar  energy  that  reaches  the  ground  in  summer  is  far  larger  than  the  net 
longwave  radiation  loss  during  the  night.  In  other  words,  much  more  heat  energy  is  supplied 
to  the  atmosphere  during  the  day  than  is  extracted  from  it  during  the  night. 
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Figiire  29.  Modeled  Horizontal  Wind  Vectors  at  14  meters  above  the  Groimd  at  (a)  0700 1st, 
(b)  14U0  1st,  and  (c)  2300  1st,  Jime  14,  1966. 


66 


The  combination  of  the  deeper  boundary  layer  and  greater  energy  available  to  heat  the 
atmosphere  during  the  day  produces  much  larger  pressure  gradients  than  are  produced  at  night. 
For  this  reason,  the  modeled  wind  directions  at  2300  1st  (Figure  29)  were  not  downslope, 
although  the  upslope  flow  components  encountered  at  1400  1st  almost  completely  disappeared. 
The  modeled  wind  direction  later  became  downslope,  but  its  magnitude  was  very  small. 

To  illustrate  diurnal  variations  of  plume  structure,  puffs  were  released  continuously  for 
29  hours  starting  at  2300  1st,  June  13,  1966.  Figure  30  shows  the  projection  of  puff  centers 
on  the  ground  at  0300,  0700,  1200,  1500,  and  2200  1st,  June  14,  1966.  The  upper-air  wind 
direction  was  southwesterly  (MI87)  for  the  entire  simulation  period,  but  puffs  remained  close 
to  the  ground  at  0300  1st  and  were  transported  to  the  southeast  The  northwesterly  flows  were 
initiated  during  the  day  when  sea-breeze  circulations  and  upslope  flows  combined. 

Nocturnal  drainage  flows  and  land  breezes  continued  to  develop  until  0700  1st  and  puffs 
were  transported  offshore.  As  the  ground  temperatures  increased  as  the  result  of  solar  heating, 
sea  breezes  and  upslope  flows  developed.  Pollutants  offshore  were  pushed  back  toward  the 
land  by  1200  1st  (Figure  30). 

Pollutants  were  initially  transported  to  the  east,  then  toward  the  south  by  the  surface  winds. 
Horizontal  wind  vectors  in  the  surface  layer  converged  into  high-altitude  locations,  where 
upward  motion  developed.  Pollutants  moved  upward  with  vertical  velocity  and  vertical  mixing 
due  to  atmospheric  turbulence.  The  pollutants  that  were  lifted  above  the  surface  boundary  layer 
encountered  southwesterly  prevailing  winds  and  moved  toward  the  northeast  as  seen  from  the 
plot  for  1500  1st. 

As  the  sun  set,  ground  temperature  continued  to  decrease  and  turbulence  subsided. 
Pollutants  remained  close  to  the  ground  and  moved  toward  the  east-southeast  at  2300  1st 
(Figure  30). 
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Figure  30.  Modeled  Trajectories  of  Puff  Centers  Projected  on  the  Siud'ace  at  (a)  0300  1st, 
(b)  0700  1st,  (c)  1200  1st,  (d)  15(K)  1st,  and  (e)  2300  1st,  J\me  14,  1966.  P\iffs 
were  released  continuously  stairting  at  2300  1st,  June  13,  1966. 
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Figuie  30.  Modeled  Trajectories  of  Puff  Centers  Projected  on  the  Surface  at  (a)  0300  1st, 
(b)  0700  1st,  (c)  1200  1st,  (d)  1500  1st,  and  (e)  2300  1st,  June  14,  1966.  Pulfs 
were  released  continuously  starting  at  2300  1st,  June  13,  1966. 
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Figure  30.  Modeled  Trajectories  of  Puff  Centers  Projected  on  the  Svirface  at  (a)  0300  1st, 
(b)  0700  1st,  (c)  1200  1st,  (d)  1500  1st,  and  (e)  2300  1st,  June  14,  1966.  Puffs 
were  released  continuously  starting  at  2300  1st,  Jvme  13,  1966. 
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Figure  30.  Modeled  Trajectories  of  Puff  Centers  Projected  on  the  Surface  at  (a)  0300  1st, 
(b)  0700  1st,  (c)  1200  1st,  (d)  1500  1st,  and  (e)  2300  1st,  June  14,  1966.  PuEs 
were  released  continuously  starting  at  2300  1st,  June  13,  1966. 
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As  seen  in  Figure  30,  surface  wind  directions  at  the  source  location  varied  360  degrees  in 
24  hours.  Accordingly,  pollutants  released  continuously  went  all  around  around  the  source. 

Figure  31  shows  the  ground  level  concentration  at  0300,  0700,  1200,  1500,  and  2200  1st, 
respectively.  The  distribution  of  the  surface  concentrations  at  0300  and  0700 1st  is  complicated, 
with  many  highs  and  lows.  These  distributions  are  definitely  not  Gaussian. 

Strong  vertical  mixing  due  to  turbulence  at  1200  and  1500 1st  resulted  in  surface  concentra¬ 
tion  distributions  similar  to  a  Gaussian  plume.  Concentration  values  decreased  monotonically 
with  the  distance  from  the  source.  Local  hot  spots  diminished  and  maximum  values  were 
much  smaller  than  the  nighttime  counterparts. 

B.  MICROVAX  2000 

The  computations  performed  by  a  Sun  4/110  workstation  were  repeated  by  a  Micro Vax 
2000  computer.  The  computational  domain  was  40  x  48  km^  with  a  horizontal  grid  spacing  of  2 
km.  Initial  and  boundary  conditions  were  identical  to  those  specified  for  the  Sun  computations. 
The  HOTMAC  simulation  started  at  0500  1st,  June  13,  1966  and  continued  for  24  hours.  In 
the  RAPTAD  simulation,  puffs  were  released  ''ontinuously  for  20  hours  stai+’*'g  at  0600  1st, 
June  13,  1966. 

We  compared  selected  HOTMAC  outputs  from  the  Sun  4/110  and  Micro  Vax  2000  simu¬ 
lations.  They  were  in  good  agreement  with  each  other.  The  differences  were,  in  gener J,  less 
than  0.1  percent.  We  consider  the  differences  insignificant 

RAPTAD  uses  the  variables  computed  by  HOTMAC  to  calculate  ptiff  transport  and  diffu¬ 
sion.  Because  a  random-number  generator  is  involved  in  RAPTAD,  instantaneous  concentration 
values  obtained  by  a  Micro  Vax  computer  are  different  from  those  obtained  by  a  Sun  worksta- 
von  even  when  HOTMAC  outputs  are  identical.  The  method  of  generating  a  random  number 
varies  from  one  computer  architecture  to  another. 
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Figxire  31.  Concentration  Distribution  (in  au'bitraxy  units)  of  Fluorescent  Peirticles  at 
Ground  Level  at  (a)  0300  1st,  (b)  0700  1st,  (c)  1200  1st,  (d)  1500  1st,  and 
(e)  2200  1st,  June  14,  1966. 
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Figure  31.  Concentration  Distribution  (in  arbitrary  units)  of  Fluorescent  Particles  at 
Ground  Level  at  (a)  0300  1st,  (b)  0700  1st,  (c)  1200  1st,  (d)  1500  1st,  and 
(e)  2200  1st,  June  14,  1966. 


Figure  32  shows  the  Micro Vax  model  results  of  puff  centers  projected  on  the  ground  at 
1100  1st,  June  13,  1966.  The  modeled  instantaneous  concentration  distributions  (in  arbittary 
units)  at  the  ground  level  are  shown  in  Figure  33. 

The  corresponding  model  results  obtained  by  a  Sun  4/1 10  workstation  are  shown  in  Figure 
34  and  Figure  35.  The  puff-center  projections  (Figure  32)  and  ground-level  concentrations 
(Figure  33)  obtained  by  a  Micro  Vax  computer  are  in  close  agreement  with  but  not  identical  to 
the  counterparts  (Figure  34  and  Figure  35)  obtained  by  a  Sun  4/1 10  workstation  because  the 
two  computers  generate  random  numbers  dilferently. 

On  the  basis  of  the  results  discussed  above,  we  conclude  that  HOTMAC  and  RAPTAD 
were  run  successfully  on  a  Micro  Vax  2000  and  that  their  results  are  in  good  agreement  with 
those  obtained  by  a  Sun  4/110  workstation. 

Finally,  we  record  in  Table  3  the  CPU  time  used  for  the  computations.  These  values 
serve  merely  as  a  reference,  as  current  workstation  technology  is  significantly  better  than  the 
computers  used  in  this  study. 


TABLE  3.  CPU  COMPARISON  BETWEEN  SUN  4/110  AND  MICROVAX  2000. 
Calculation  Sun  4/1 10  Micro  Vax  2000 


HOTMAC 

21  X  25  X  16  grid  points  4  h,  11  min  22  h,  10  min 

28-h  simulation 

RAPTAD  26  min  3  h,  45  min 

20-h  simulation 
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Figure  32.  Modeled  Trajectories  of  Puff-Centers  Projected  on  the  Surface  at  1100  1st, 
Jui  e  13,  iOGG.  Puffs  were  released  continuously  for  20  hours  starting  at 
0600  kt,  June  13,  1966.  The  resxilts  were  obtained  by  a  MicroVax  2000. 
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Figure  34.  Same  as  in  Figure  32,  but  the  Results  were  Obtained  by  a  Stm  4/110  Work¬ 
station. 
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SECTION  vn 


BUOYANT  PLUMES 

In  this  section,  we  will  investigate  positive  and  negative  buoyancy  effects  in  HOTMAC 
and  RAPTAD  (Phase  IH  tasks). 

The  most  widely  used  method  in  treating  a  buoyant  plume  is  based  on  a  plume-rise  theory 
advanced  by  Briggs  (Reference  14).  One  can  estimate  the  height  where  the  plume  levels  off. 
This  height  is  referred  to  as  an  effective  stack  height  Therefore,  the  simplest  way  to  deal  with 
a  buoyant  plume  is  to  apply  the  Briggs  formula  in  the  initial  stage  when  the  buoyancy  effect 
is  significant.  When  the  plume  becomes  neutrally  buoyant,  RAPTAD  can  be  applied  as  before 
except  that  a  stack  height  is  replaced  by  an  effective  stack  height 

There  have  been  some  efforts  (e.g..  Reference  15)  to  incorporate  the  buoyancy  effect  into 
a  puff  model.  Thus,  a  puff  model  can  be  used  firom  an  initial  stage  and  a  patchwork  with  a 
plume-rise  model  is  not  needed.  Indeed,  the  puff  model  by  Gaffen  et  al.  (Reference  15)  used 
a  concept  very  similar  to  that  used  for  a  plume-rise  model  to  include  the  buoyancy  effect 

The  advantage  of  a  plume-rise  model  is  its  simplicity.  A  plume  rises  because  its  density 
is  smaller  than  that  of  the  ambient  air,  and  it  levels  off  when  it  loses  buoyancy  by  mixing 
with  the  ambient  air.  The  entrairunent  process,  ot  mixing  rate,  is,  however,  dependent  upon 
a  complex,  turbulent  heat  and  mass  exchange  between  a  puff  and  the  ambient  air,  but  the 
entrainment  process  is  greatly  simplified  in  a  plume-rise  model.  Both  a  plume-rise  model  and 
a  puff  model  assume  that  a  buoyant  plume  does  not  perturb  the  dynamics  of  the  ambient  flows. 
In  other  words,  the  ambient  flow  conditions  influence  the  plume  dynamics,  but  not  vice  versa. 
For  this  reason,  a  plume-rise  model  and  a  puff  model  are  referred  to  as  a  “passive  plume-rise 
model”  in  this  study. 
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On  the  other  hand,  we  expect  that  a  high-intensity  buoyant  plume,  such  as  one  that  would 
result  from  a  large  fire  or  explosion,  would  influence  significantly  the  dynamics  of  the  ambient 
flow  conditions.  A  fire  heats  the  air  surrounding  it  and  produces  a  low-pressure  area.  Horizontal 
pressure  gradients  generated  by  a  fire  force  the  air  to  move  toward  the  fire.  Convergence  of 
horizontal  wind  vectors  over  a  fire  creates  a  vertical  motion  or  a  plume  rise.  Therefore,  wind 
and  temperature  distributions  are  significantly  different  with  and  without  a  large  heat  source. 

Furthermore,  a  large  fire  generates  vigorous  turbulent  motion  because  temperature  strati¬ 
fication  becomes  extremely  unstable  as  a  result  of  intense  heating. 

A  passive  plume-rise  model  cannot,  and  is  not  designed  to,  address  the  alterations  of  the 
ambient  flow  conditions.  Thus,  its  application  is  limited  to  a  relatively  low  intensity  heat  source. 

An  equally  interesting  problem  arises  when  a  high-density  and/or  a  low-temperature  gas 
is  released  in  the  air.  If  a  dense  gas  is  spilled  at  the  surface,  the  gas  spreads  horizontally  and 
can  stay  in  the  surface  layer  for  a  long  period  of  time,  because  the  density  stratification  over  a 
contaminated  area  is  very  stable  and  turbulent  mixing  is  substantially  suppressed.  Therefore, 
we  expect  that  a  dense  gas  behaves  and  impacts  the  ambient  flows  quite  differently  from  a 
buoyant  plume.  For  this  reason,  people  have  developed  different  models  to  treat  a  dense  gas 
and  a  plume  rise. 

In  this  section,  we  will  investigate  the  feasibility  of  incorporating  positive  and  negative 
source  buoyancy  effects  in  HOTMAC  and  RAPTAD.  Our  approach  is  quite  different  from  that 
of  a  passive  plume-rise  model.  We  specify  a  heat  source  as  a  boundary  condition;  then  the 
governing  equations  in  HOTMAC  compute  wind,  temperature,  and  turbulence  distributions. 
RAPTAD  uses  the  HOTMAC  variables  to  determine  the  centroid  location  and  size  of  each 
puff.  Unlike  a  plume-rise  model,  a  puff  in  RAPTAD  is  a  tracer  that  follows  precisely  the 
turbulent  fluid  motion.  Our  approach  in  this  study  is  referred  to  as  a  “dynamic  plume  model” 


86 


The  advantage  of  a  dynamic  plume  model  is  that  the  fluid  motions  and  turbulent  processes 
are  consistent  with  the  physical  laws  expressed  by  the  governing  equations.  In  other  words, 
we  can  minimize  the  ambiguities  and  inconsistencies  introduced  in  a  passive  plume  model. 
For  example,  methods  used  in  a  passive  plume  model  to  define  the  parameters  of  turbulent 
entrainment  processes  are  based  on  simple  and  idealized  conditions.  A  passive  plume  model, 
by  definition,  cannot  provide  feedback  to  the  ambient  flow  conditions.  For  example,  vertical 
motions  resulting  from  a  plume  rise  do  not  affect  the  distribution  of  horizontal  wind  components 
in  the  background  flow.  Of  course,  this  is  inconsistent  with  mass  continuity. 

The  disadvantage  of  a  dynamic  plume  model  is  its  complexity.  We  must  solve,  numer¬ 
ically,  a  set  of  three-dimensional  hydrodynamic  models  with  appropriate  initial  and  boundary 
conditions.  Reaching  these  solutions  requires  considerable  knowledge  and  experience  in  fluid 
dynamics  and  numerical  modeling.  Furthermore,  the  computations  could  be  quite  expensive, 
because  a  very  small  integration  time  step  is  used  to  ensure  computational  stability. 

Therefore,  a  dynamic  plume  model  is  currently  limited  to  research  applications.  We 
expect  that  it  will  become  practical  for  an  operational  use  in  the  near  future,  because  computer 
capabilities  relative  to  cost  will  increase  substantially  in  the  future,  particularly  for  engineering 
workstations. 

We  now  discuss  results  of  simulations  from  HOTMAC  and  RAPTAD  where  positive  and 
negative  buoyant  sources  are  placed  at  the  lower  boundary  of  a  computational  domain.  We 
use  an  area  heat  source  to  simulate  plume  rise  associated  with  a  fire  and  an  area  cold  air 
distribution  to  simulate  dense  gas  dispersion  in  the  surface  layer.  Results  are  preliminary, 
because  our  purpose  is  to  test  the  feasibility  of  incorporating  a  large  heat  source  or  a  sink  into 
HOTMAC.  We  did  not  attempt  to  evaluate  the  model  performance  against  any  observations 
in  this  study. 
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A.  POSITIVE  BUOYANT  SOURCE 


We  have  selected  the  simplest  boundary  condition  because  we  would  like  to  isolate  the 
effects  of  the  positive  buoyancy  source  on  air  flow  from  any  other  causes  that  might  influence 
the  flow.  We  have  chosen  a  flat  surface  boundary  and  a  nighttime  simulation  period  to  exclude 
the  diurnal  variation  of  the  solar  heating. 

The  computational  domain  was  20  x  20  x  10  km  in  x,  y,  and  z  directions,  respectively, 
and  the  horizontal  grid  spacing  was  500  meters.  The  grid  spacing  in  the  vertical  direction  is 
variable  with  height:  4  meters  near  the  surface  and  1000  meters  close  to  the  top  of  computational 
domain.  A  total  of41x41x31  grid  points  represented  the  computational  domain.  An  area 
heat  source  of  3  x  4  km^  was  placed  between  x  =  14.25  and  17.25  km  and  y  =  7.75  and  1 1.75 
km  in  the  computational  domain.  We  adopted  a  heat  flux  of  50,000  W/m^  ,  which  corresponds 
to  an  extremely  large  fire. 

Initially,  vrind  and  potential  temperature  distributions  were  assumed  uniform  in  the 
horizontal  directions.  The  vertical  profile  of  the  irutial  east-west  wind  component  was  obtained 
from  a  logarithmic  profile  U/u.  =  1/k  In  (z  +  Zo)/zo  where  u.  =  0.2  m/s,  Zo  =  0.1  meters,  and 
k  =  0.4.  The  north-south  wind  component  was  zero  iiutially.  The  vertical  profile  of  the  irutial 
potential  temperature  was  20°C  at  the  ground  and  increased  with  height  at  a  rate  of  0.001  °C/m 
up  to  3000  meters  above  msl  and  at  a  rate  of  0.003®Cym  above  that  level. 

Wind  fields  appear  to  reach  a  quasi-steady  state  at  20  minutes  after  the  initiation  of  heating. 
The  heating  rate  was  kept  constant  during  the  entire  simulation  period.  Thus,  wind,  temperature, 
and  turbulence  distributions  discussed  here  are  for  those  at  20  minutes  after  the  fire  was  initiated. 
Figure  36  shows  a  horizontal  wind  vector  distribution  at  30  meters  above  the  ground.  Horizontal 
wind  vectors  converge  over  a  point  slightly  downwind  of  the  western  boundary  of  the  fire.  If 
the  prevailing  wind  were  zero,  the  convergence  should  have  occtured  at  the  center  of  the  fire. 
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Figure  36.  Horizontal  Wind  Vector  Distribution  at  30  meters  above  the  Ground.  The 
heated  area  is  indicated  by  a  rectangle. 
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The  convergence  of  the  horizontal  wind  vectors  decreased  in  magnitude  with  height,  and  the 
wind  vector  distribution  switched  to  divergence  at  the  higher  levels.  This  is  the  consequence 
of  the  mass  continuity  constraint  imposed  in  the  model.  The  switch  from  convergence  to 
divergence  occurred  at  approximately  2000  meters  above  the  ground.  Perturbations  caused  by 
a  fire  were  negligible  at  levels  higher  than  6500  meters  above  the  ground. 

Wind  vectors  in  east-west  vertical  cross  sections  along  the  southern  boundary  and  the 
central  axis  of  the  fire  are  shown  in  Figure  37.  These  wind  vectors  show  large  undulations: 
a  strong  subsidence  in  front  of  the  fire  and  an  extremely  large  lift  over  the  convergence  zone 
followed  by  another  subsidence.  Convergence  of  horizontal  wind  vectors  in  the  surface  layer 
and  retiun  flows  generated  weU-estabUshed  rotors  downwind  of  the  fire.  The  maximum  vertical 
velocity  was  43.5  m/s,  and  the  maximum  horizontal  velocity  was  18.5  m/s. 

Wind  vectors  in  north-south  vertical  cross  sections  along  x  =  14.5  km,  16  km,  and  17 
km  are  shown  in  Figure  38.  Wind  vector  distributions  in  the  vertical  cross  section  along  the 
western  boundary  of  the  fire  (x  =  14.5  km)  clearly  show  almost  symmetric  vortices  for  each 
side  of  the  fire.  The  wind  vectors  along  the  north-south  centerline  of  the  heated  area  (x  =  16 
km)  exhibit  much  smaller  vortices  than  those  at  x  =  14.5  km.  The  updraft  over  the  central 
part  of  the  heated  area  diminishes  at  approximately  3.5  km  above  the  ground  and  is  replaced 
by  a  downdraft  in  the  upper  levels.  Wind  vectors  along  the  eastern  boundary  of  the  fire  (x  = 
17  km)  indicate  that  the  air  motions  are  mostly  downward. 

The  fire  produced  very  large  vertical  velocities  (over  40  ra/s)  along  the  fire  column  over  the 
convergence  point  of  the  horizontal  wind  vectors.  The  question  arises  whether  the  hydrostatic 
approximation  assumed  here  is  stil^  valid.  We  repeated  the  simulation  with  a  nonhydrostatic 
code  and  obtained  results  very  similar  to  those  reported  here.  The  magnitude  of  the  dynamic 
pressure  obtained  in  the  nonhydrostatic  pressure  computation  was  less  than  1  percent  of  the 
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Figure  37.  Wind  Vectors  in  East- West  Vertical  Cross  Sections  along  (a)  the  Southern 
Boundary  and  (b)  the  Central  Axis  of  the  Fire. 
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Figtire  37.  Wind  Vectors  in  East-West  Vertical  Cross  Sections  along  (a)  the  Southern 
Boundary  and  (b)  the  Central  Axis  of  the  Fire. 
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Figure  38.  Wind  Vectors  in  the  North-South  Vertical  Cross  Sections  along  (a)  x  =  14.5  km, 
(b)  X  =  16  km,  and  (c)  x  =  17  km. 
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Figure  38.  Wind  Vectors  in  the  North-South  Vertical  Cross  Sections  along  (a)  x  =  14.5  km, 
(b)  X  =  16  km,  and  (c)  x  =  17  km. 
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Figure  38.  Wind  Vectors  in  the  North-South  Vertical  Cross  Sections  along  (a)  x  =  14.5  km, 
(b)  X  =  16  km,  and  (c)  x  =  17  km. 
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hydrostatic  pressure.  Thus,  we  tentatively  conclude  that  the  hydrostatic  pressure  approximation 
is  still  valid  even  for  a  large  heat  source.  This  conclusion  may  depend  on  the  horizontal  grid 
resolution  used  in  a  model  because  a  coarse  grid  tends  to  underestimate  horizontal  derivatives 
of  wind  distribution.  Those  terms  are  the  main  contributors  to  the  dynamic  pressure. 

The  fire  modified  considerably  not  only  mean  winds  but  also  turbulence  variables.  A  heat 
flux  of  50,000  W/m^  over  the  fire  area  resulted  in  a  maximum  surface  temperature  of  367°  C 
near  the  western  boundary  of  the  fire  where  the  convergence  of  horizontal  wind  vectors  was 
observed.  The  surface  temperatures  over  the  heated  area  were  over  300°  C  everywhere.  The 
heat  energy  from  a  fire  becomes  a  buoyancy  source  of  atmospheric  turbulence.  Figure  39 
shows  the  vertical  profiles  of  standard  deviation  of  vertical  velocity,  twice  turbulence  kinetic 
energy,  and  standard  deviation  of  potential  temperature  at  site  1  (Figure  36),  which  is  located 
near  the  center  of  the  western  boundary  of  the  fire.  The  values  for  turbulence  kinetic  energy 
reported  in  Figure  39  are  at  least  20  times  larger  than  those  encoimtered  in  the  convective 
atmospheric  boundary  layer.  Turbulence  reached  6000  meters  above  the  ground,  whereas  a 
typical  atmospheric  boundary  layer  height  is  approximately  2000  meters. 

To  visualize  the  extent  of  a  plume  rise,  we  released  puffs  continuously  over  the  heated 
area.  Figure  40  shows  three-dimensional  projections  of  puff  centroids  at  20, 40,  and  60  minutos 
after  the  release,  respectively.  Puffs  close  to  the  ground  first  moved  toward  the  convergence 
point,  which  is  located  near  the  center  of  the  western  boundary  of  the  fire.  The  updraft  over 
the  convergence  zone  transported  puffs  upward  over  8800  meters  above  the  ground.  Puffs 
detrained  from  the  fire  coluiim  in  the  levels  between  2000  and  88(X)  meters  above  the  ground 
where  horizontal  wind  vectors  were  divergent 
B.  NEGATIVE  BUOYANT  SOURCE 

Because  of  the  nature  of  the  chemicals  dealt  with  at  VAFB,  little  emphasis  is  placed 
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Figure  39.  Vertical  Profiles  of  Standard  Deviation  of  Vertical  Velocity,  Twice  Turbu¬ 
lence  Kinetic  Energy,  q^;  and  Standard  Deviation  of  Potential  Temperatvire, 
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Figure  40.  Three-Dimensional  Projections  of  PulF  Centroids  at  (a)  20  minutes,  (b)  40  min¬ 
utes,  and  (c)  60  minutes  after  the  Release. 
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Figure  40.  Three-Dimensional  Projections  of  Puff  Centroids  at  (a)  20  minutes,  (b)  40  min 
utes,  and  (c)  60  minutes  after  the  Release. 


on  a  negative  buoyant  source  at  this  initial  stage  of  testing  the  feasibility  of  operating 
HOTMAC/RAPTAD  with  the  VAFB  terrain.  Gas  can  be  heavier  than  air  when  its  molecular 
weight  is  greater  than  that  of  air,  or  when  its  temperature  is  lower  than  that  of  air,  or  when 
both  conditions  exist.  In  this  section,  the  terms  “negative  buoyant  source,”  “gas  heavier  than 
air,”  and  “dense  gas”  are  used  interchangeably. 

The  approach  used  to  incorporate  a  negative  buoyant  source  in  HOTMAC  is  similar  to 
that  used  to  incorporate  a  positive  buoyant  source,  as  discussed  in  the  previous  subsection. 
Again,  we  made  no  effort  to  verify  the  model  results  with  observations. 

We  selected  a  computational  domain  of  80  x  140  x  5  km  in  x,  y,  and  z  directions, 
respectively.  For  simplicity,  we  assumed  that  the  surface  was  flat  and  that  wind  and  temperature 
distributions  were  initially  uniform  in  the  horizontal  directions.  Initially,  wind  direction  was 
southerly  and  wind  speed  was  3  m/s.  Initial  potential  temperature  was  30®C  at  the  surface  and 
increased  with  height  at  a  rate  of  0.010°C/m  up  to  700  meters  above  the  ground.  The  lapse  rate 
decreased  to  0.00236'’C/m  for  the  rest  of  the  levels  up  to  the  top  of  the  computational  domain. 

Horizontal  grid  spacing  was  4  km,  and  vertical  grid  spacing  varied  with  height.  We  placed 
near  the  center  of  the  computational  domain  a  cold  area  of  8  x  8  km  with  a  temperature  of 
-150°C.  The  integration  was  initiated  at  0300  1st  and  lasted  for  6  hours. 

Figure  41  shows  the  distribution  of  horizontal  wind  vectors  at  10  meters  above  the  ground 
at  0400 1st.  As  expected,  the  wind  field  diverges  over  the  cold  area  where  the  surface  pressure 
was  higher  than  the  surroundings.  The  wind  perturbations  were  confined  in  the  layer  less 
than  30  meters  ffan  the  ground,  and  turbulence  intensity  was  very  small.  We  released  puffs 
continuously,  starting  at  0400  1st  Figure  42  shows  the  locations  of  puff  centroids  projected  on 
the  grotmd  at  0500  1st  and  0900  1st,  respectively.  Puffs  were  initially  coherent  but  bifurcated 
with  time  because  of  the  divergent  nature  of  the  wind  field  over  the  cold  area  (Figure  41). 


101 


(  1  t  1  t  t  1  1  1  t  t  t  t  t  1 

t  t  \  t  t  t  t  t  1  t  I  t  t  t  1 

t  \  t  t  t  (  I  t  1  t  t  t  t  t  « 

\  t  t  t  t  t  t  t  t  t  t  1  t  I  t 

\  t  t  t  t  t  t  t  t  t  t  I  I  I  t 

t  \  \  t  1  1  t  t  1  t  t  t  I  I  t 

t  \  \  1  1  t  t  t  t  t  t  t  }  f  t 

\  \  \  \  1  i  1  t  t  I  t  I  f  I  t 

N\\\N\\tt////// 

\NVS\\\  \  1  t  t  f  /  ^  * 

\  S  V  N  N  N  \  X  >  »  flB’  ^  ' 

\NVVNNK  \  %  #  ^  ^  ^  ^  4 

^  m  ^  4 

^  ^  ^  ^  ^  ^  9  ^  %  «  »  •  »  # 

I  •  «  •  *  • 

•  «  «  «  t 

•  «  •  %  I  t 


»  ^  # 


'  •  .  0]r 


4  0%% 


\ 

\ 

\ 

\ 

\ 

% 

% 

s 

% 

« 

t 

% 

\ 

\ 

\ 

\ 

\ 

\ 

\ 

\ 

% 

« 

\ 

\ 

\ 

> 

\ 

( 

\ 

\ 

1 

X 

x 

X 

\ 

\ 

1 

X 

X 

X 

X 

\ 

X 

X 

X 

X 

\ 

\ 

X 

X 

X 

X 

X 

X 

X 

X 

X 

\ 

\ 

X 

X 

X 

\ 

t 

X 

X 

X 

X 

X 

\ 

\ 

\ 

X 

X 

\ 

\ 

X 

X 

x 

) 

X 

\ 

x 

X 

X 

X 

X 

X 

X 

X 

X 

t 

X 

\ 

x 

x 

X 

X 

X 

\ 

X 

X 

X 

X 

X 

\ 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

Figure  41.  Horizontal  Wind  Vectors  at  10  meters  above  the  Ground  at  0400  1st.  The  cold 
area  is  indicated  by  a  rectangle. 
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Figure  42.  Projections  of  Puff  Centroids  on  the  Ground  at  (a)  0500  1st  and  (b)  0900  1st. 
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Figure  42.  Projections  of  Puff  Centroids  on  the  Ground  at  (a)  0500  1st  and  (b)  0900  1st. 
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Figure  43  shows  the  surface  concentration  contours  (in  arbitrary  units)  at  0500 1st  and  09(X) 
1st,  respectively.  Both  figures  show  that  the  distributions  are  not  Gaussian,  and  bifurcations 
are  evident,  particularly  at  0900  1st. 
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Figure  43.  Surface  Concentration  Contours  (in  arbitrary  units)  at  (a)  0500 1st  and  (b)  0900 1st. 
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Figure  43.  Surface  Concentration  Contours  (in  arbitrary  units)  at  (a)  0500 1st  and  (b)  0900 1st. 


SECTION  vm 


CONCLUSIONS  AND  RECOMMENDATIONS 


A.  CONCLUSIONS 

•  It  is  feasible  to  operate  the  three-dimensional  atmospheric  models  HOTMAC  and 
RAPTAD  to  forecast  the  transport  and  dispersion  of  airborne  materials  at  VAFB. 

•  HOTMAC  and  RAPTAD  predictions  were  at  least  as  good  as  those  obtained  by  diagnostic 
models  and  the  MI  model  where  wind  data  were  available. 

•  HOTMAC  and  RAPTAD  predictions  were  far  better  than  those  obtained  by  diagnostic 
models  and  they  were  the  best  practical  solution  where  wind  data  were  not  available. 

•  It  is  feasible  to  operate  HOTMAC  and  RAPTAD  on  a  desktop  computer.  HOTMAC  took 
4  hours  11  minutes  and  22  hours  10  minutes  CPU  time,  respectively,  on  a  Sun  4/100 
workstation,  and  a  MicroVax  2000  computer  for  a  28  hour  forecast  with  21  x  25  x  16 
grid  points.  RAPTAD  took  26  minutes  and  3  hours  45  minutes  CPU  time,  respectively, 
on  a  Sun  4/1 10  workstation,  and  a  MicroVax  2000  computer  for  a  20  hour  simulation. 

•  The  affordability  and  portability  of  the  desktop  computer  has  opened  the  door  to 
upgrading  toxic  hazard  modeling  capabilities  for  emergency  response  management  at 
VAFB. 

•  The  HOTMAC  and  RAPTAD  modeling  system  would  be  a  useful  addition  to  enhance 
current  VAFB  emergency  response  management  capabilities. 

•  It  is  feasible  to  incorporate  the  positive  and  negative  source  buoyancy  effects  in  HOT¬ 
MAC  and  RAPTAD.  Our  approach  can  minimize  the  ambiguities  and  inconsistencies 
introduced  in  a  simple  plume  rise  model. 
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B.  RECOMMENDATIONS 


•  Become  familiar  with  background  theories  and  operating  procedures  of  HOTMAC  and 
RAPTAD.  Install  the  models  on  a  mainframe  computer  and  repeat  the  simulations 
discussed  in  this  report 

•  Consider  upgrading  computer  capabilities  for  emergency  response  applications  at  VAFB. 
MicroVax  computers  may  not  be  fast  enough  to  operate  HOTMAC  and  RAPTAD. 

•  Modify  HOTMAC  to  run  all  the  time  and  store  wind  and  turbulence  data  for  the  next  24 
hours.  When  an  emergency  occurs,  RAPTAD  can  be  used  immediately  with  the  wind 
data  stored  on  a  disk.  The  RAPTAD  computation  is  much  faster  than  that  of  HOTMAC. 
Thus,  this  approach  meets  better  the  time  constraint  imposed  under  emergency  situations. 

•  Develop  a  method  to  integrate  tower  data  into  HOTMAC.  The  wind  data  can  be  used  to 
initialize  and  correct  the  wind  distribution  in  HOTMAC.  These  may  be  accomplished  by 
a  dynamical  initialization  technique  and  a  four-dimensional  data  assimilation  method. 

•  Develop  a  method  to  predict  conce:  ration  variances  which  can  be  used  to  estimate 
uncertainties  associated  with  predictions.  A  short  time  averaging  value  is  required 
for  predicting  concentration  of  toxic  materials.  Such  a  value  normally  exhibits  great 
variations  in  time  and  space. 

•  Add  model  physics  necessary  to  simulate  the  evolution  of  fog  formation  and  dissipation 
processes.  Fog  is  frequently  observed  at  VAFB  and  it  affects  the  heat  energy  balance 
at  the  grour  * 

•  Continue  to  investigate  the  feasibility  of  incorptHating  positive  and  negative  source 
buoyancy  effects  in  HOTMAC  and  RAPTAD  and  test  the  scheme  with  observations. 

(The  reverse  of  this  page  is  blank) 
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APPENDIX  A 


DESCRIPTION  OF  MODELS  AND  NUMERICAL  PROCEDURES 
A.  HOTMAC  (Higher  Order  Tiirbulence  Model  for  Atmospheric  Oirciilation) 

1.  Model  Equations 

This  model,  also  referred  to  as  a  “second-moment  turbulence-closure  model,”  is 
beised  on  a  set  of  second-moment  tvubulence  equations  closed  by  assuming  certain  rela¬ 
tionships  between  unknown  higher-order  turbulence  moments  and  the  known  lower-order 
variables.  The  model  output  variables  are  winds,  potential  temi>erature,  mixing  ratios 
of  water  vapor  and  liquid  water,  turbulence  second  moments,  a  tvubulence  length  scale, 
and  turbulence  transport  coefficients  (eddy  viscosity  and  eddy  diffusivity).  These  results 
can  be  used  as  inputs  to  pollutant  dispersion  models.  The  model  is  time  dependent  and 
three-dimensional  in  space. 

HOTMAC  can  be  used  under  quite  general  conditions  of  flow  and  thermal  strati¬ 
fication;  methods  for  tiubvdence  parameterization  are  more  advanced  than  those  in  simple 
eddy  viscosity  models.  The  present  model,  combined  with  a  statisticed  cloud  model,  has 
simvilated  interaction  between  water  phase  changes  and  basic  dynamic  variables.  For  exeun- 
ple,  computed  tvubulence  energy  increases  substantially  in  the  layers  where  condensation 
occvus.  This  appears  reasonable,  since  the  latent  heat  released  by  condensation  produces 
loc2il  vmstable  layers,  resulting  in  generation  of  tvubulence.  Effects  of  short-  and  long-wave 
solar  radiation,  tall  tree  canopies,  and  topography  are  also  included  in  the  model.  Svuface 
temperatvues  are  computed  &om  a  heat  conduction  equation  for  the  soil  and  a  heat  energy 
balance  equation  at  the  svuface. 

The  present  model  assumes  hydrostatic  equilibrimn  and  uses  the  Boussinesq 
approximation.  Therefore,  in  theory,  the  model  applications  are  limited  to  flows  where  the 
loczd  acceleration  emd  advection  terms  in  the  equation  of  vertical  motion  are  much  smaller 
than  the  acceleration  due  to  gravity  (hydrostatic  equilibrium)  and  temperatvue  variations 
in  the  horizontal  are  not  too  large  (Boussin^q  approximation). 
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The  model  has  been  used  for  a  variety  of  fluid  problems:  surface  boimdary  layer 
(Reference  16),  atmospheric  boundary  layer,  (Reference  17;  Reference  18;  Reference  19; 
Reference  20;  Reference  21),  airflow  over  tall  tree  canopies  (Reference  22),  air  pollution 
transport  (Reference  23),  ocean  boimdary  layer,  (Reference  24),  airflow  over  a  cooling  pond 
(Reference  25),  laboratory  flows  (Reference  26),  flow  over  complex  terrain  (Reference  27), 
and  the  results  have  been  used  in  dispersion  simulations  over  complex  terredn  (Reference  28; 
Reference  29;  Reference  8),  over  the  eastern  half  of  the  U.S.  (Reference  30)  and  over 
the  mountainous  western  U.S.  (Reference  31).  This  model  has  also  been  used  by  others 
including  Dobosy  (Reference  32),  Shaw  (Reference  33),  Burk  (Reference  34),  Miyakoda 
and  Sirutis  (Reference  35  ),  Freeman  (Reference  36),  and  Sim  and  Ogura  (Reference  37). 
A  complete  summary  of  the  model  is  given  in  a  recent  review  paper  (Reference  38). 

The  basic  equations  of  HOTMAC  for  mean  wind,  temperature,  mixing  ratio  of 
water  vapor,  and  turbulence  are  similar  to  those  used  by  Yamada  in  Reference  28  and 
Reference  39,  except  two  improvements,  the  nested  grid  capability  and  effects  of  shadows 
produced  by  terrain,  are  added. 

A  terrain-following  vertical  coordinate  system  is  used  in  order  to  increase  the 
accuracy  in  the  treatment  of  surface  boimdary  conditions: 

=  ,  (^-1) 

where  z*  and  z  are  the  tremsformed  and  C8u1;esi2m  vertical  coordinates,  respectively;  Zg  is 
ground  elevation;  H  is  the  material  surface  top  of  the  model  in  the  z*  coordinate;  and  H 
is  the  corresponding  height  in  the  z-coordinate. 

The  governing  equations,  following  the  coordinate  transformation,  are  (Refer¬ 
ence  28) 
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5r  V  ^  dx  J  5t/  \  dy  )  ^  H  —  Zg  dz* 


(—uw) 


{A -2) 


112 


DV 

Di 


<  >\  dZg 

Qv  )  dy 


+ 


H 

H-Zg 


d 

dz* 


(—vw) 


{A -3) 


_ }^(u^  +  v^)  =0 

dx  dy  dz*  H  —  Zg\  dx  dy ) 

where 


and 


W*  = 


H 

H-Zg 
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z*-H 

H-Zg 
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Di  -  dt  ^  dx  ^  dy  ^  dz* 


(A -6) 


In  the  above  expressions,  <>  indicates  an  average  over  a  horizontal  surface.  The  second 
terms  on  the  right-hand  side  of  Equations  (A-2)  and  (A-3)  indicate  the  effects  of  groimd 
slope.  For  simplicity,  H  is  specified  as 


H  =  H  +  z 


gmax  j 


(A -7) 


where  Zgmax  is  the  maximtun  value  of  the  ground  elevation  in  the  computational  domain. 
The  geostrophic  winds  Ug  and  Vg  are  computed  from  Reference  28,  i.e.. 


fU,  =  fV,(H) 


H-Zg  1  d  ^  , 

<Q^{H)>^^  H  <er>>dy^  ”  ^ 

9  dzg  >  / 

H  dy  J„  <  > 


(A -8) 


and 
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g  dzg  f  AGf,  , 

H  dx  J,,  <Q,>  "  ’ 


— 1 — 
<Qv  >  dx 


{A -9) 


where  A0„  =  0„—  <  0„  >,  and  the  abbreviated  symbols  Ug{H)  =  Ug{x,y,H,t). 
Vg{H)  =  Vg{x,y,H,t)^  axe  used.  The  derivations  of  Elquations  (A-2)  to  (A-5),  (A-8),  and 
(A-9)  axe  given  in  the  Appendix  of  Reference  28.  A  turbulence  kinetic  energy  equation  is 
given  by 

Dt\2  )  ax  [  *5xV2yJ  ay  [  V2yJ 
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+  0gwdv 
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2tnd  a  turbulence  length  scale  i  is  obtauned  from 


(A  - 10) 


-  [ft;; 


(A -11) 


where  +  v^  +  is  twice  the  turbulence  kinetic  energy,  w9^  turbulence  heat  flux, 

$v  the  fluctuation  part  of  virtual  potential  temperature,  and  {Fi,F2,Sq,Si,  and  Bi)  = 
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(1.8,  1.33,  0.2,  0.2,  eind  16.6),  empiriceil  constants  determined  from  laboratory  experiments 
(Reference  38).  The  internal  heat  energy  equation  is  written  as  in  Reference  8 
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(A -12) 


The  long- wave  radiation  flux  Rn ! pCp  is  computed  according  to  Sasamori  (Reference  40). 
A  conservation  equation  for  mixing  ratio  of  water  vapor  is  given  by 


DQv 
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dQv-i 
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H  —  Zg  dz* 


(-wqv) 


(A  -  13) 


The  turbtdent  fluxes  in  Equations  (A-2),  (A-3),  (A-10),  (A-11),  (A-12),  and  (A-13)  axe 
obtained  from  simplified  second-moment  turbulence-closure  equations  (Reference  27): 


{uw,  vw)  =  —£qS M 
{w9,w^)  =  -aiq^M 


dz  ’  dz 

dB  dQv 
dz'  dz 


{A  —  14a,  b) 
{A  —  15a,  b) 


where  ^  m  and  oc  are  fimctions  of  the  flux  Richardson  number,  and  q(=  KhJKm  where  Kh 
is  an  eddy  diffusivity  coefficient  and  Km  is  an  eddy  viscosity  coefficient)  is  the  reciprocal 
of  the  turbulent  Prandtl  number. 

The  expressions  for  'S  m  and  a  were  obtaiined  from  the  Level  2  model  of  Mellor 
£md  Yamada  (Reference  17)  where  temporal  and  spatizd  derivatives  in  Equation  (A-10)  are 
neglected.  The  readers  are  referred  to  Reference  19  for  further  discussions  of  the  Level  2 
model  equations.  The  final  expressions  for  ^  m  and  a  are  given  in  Reference  27  and  are 
not  repeated  here. 
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2.  Boundaxj"  Conditions 

Surface  boundary  conditions  for  Equations  (A-2),  (A-3),  and  (A-10)  to  (A-13) 
are  constructed  from  the  empiricaJ  formulas  by  Dyer  and  Hicks  (Reference  10)  for  the 
nondimensional  wind  emd  temperature  profiles  (see  the  Appendix  of  Reference  28).  Strictly 
speaking,  the  formulas  are  valid  only  for  horizontally  homogeneous  surfaces.  It  is  assumed, 
however,  that  the  same  relations  axe  fair  approximations  over  nonhomogeneous  terrain, 
provided  that  the  formulas  are  apphed  sufficiently  close  to  the  surface.  It  should  be  noted 
that  vegetation  plays  an  active  part  in  the  apportionment  of  available  heat  energy  between 
convective  (sensible  and  latent)  and  conductive  (into  the  soil)  components.  Use  of  the 
similarity  formulas  requires  knowledge  of  the  surface  temperattires;  a  method  to  obtain 
the  surface  temperature  is  discussed  below. 

The  temperature  T,  in  the  soil  layer  is  obtained  by  solving  the  heat  conduction 

equation 


dt  ~dz,\  ^dzj  ’ 


(A  - 16) 


where  z,  is  positive  downward,  and  soil  difFusivity  K,  can  be  a  function  of  soil  moisture 
content.  Appropriate  bovmdary  conditions  for  solution  of  Equation  (A- 16)  are  the  heat 
energj'  balance  at  the  soil  surface  and  specification  of  the  soil  temperature  or  soil  heat  flux 
at  a  certain  depth  whose  value  is  dependent  on  the  duration  of  the  integration.  The  heat 
energy  badance  at  the  surface  is  given  by 


Rs  +  Rl  i-RLl=H,  +  LE  +  Gs  ,  (A -17) 

where  R,  is  the  incoming  direct  solar  radiation  absorbed  by  the  surface,  J.  is  the 
incoming  long- wave  radiation,  and  Ri is  the  outgoing  long- wave  radiation. 

The  surface  heat  flux  H,,  latent  heat  flux  LE,  and  groxmd  heat  flux  G,  are  given 
by 


R$  —  j 


(A- 18) 
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LE  —  , 


{A -19) 


and 

,  {A  -  20) 

G 

where  pa  is  the  air  density,  tt*  is  the  friction  velocity,  T*  is  the  temperature  scale, 
is  the  water  vapor  scale,  and  the  subscript  G  denotes  the  value  at  the  groimd  surface. 
Substituting  Equations  (A-18)  to  (A-20)  into  Equation  (A-17)  we  obtain 


az. 


Ra  + eRi  i -eaTc  -PaCpuJT^ilA  B  ^)  -  KadTsldza\G  ,  (A -21) 


where  the  relation 


Rl  T=  +  (1  -  t)Ri  i  ,  (A  -  22) 

and  Bowen  ratio 


B  =  HJLE  ,  (A  -  23) 

are  used;  e  is  the  emissivity  of  the  sxirface  and  <7  is  the  Stefan-Boltzman  constant.  Garratt 
and  Hicks  (Reference  41)  obtained  a  relationship  between  the  STirface  temperature  and  air 
temperature  at  zi  (in  the  surfaw:e  layer); 


(€>(ri)  -  ©g)T,  ^  =  (Pr/A:)[(£n{(2i  +  Zot)lzo}  +  tn{zo!zot)  -  V’]  ,  (A  -  24) 

where  Pr  is  the  turbvdence  Prandtl  number  at  neutral  stability,  k  is  the  von  Karman 
constant,  2<,  and  Zot  are  the  roughness  lengths  for  momentum  and  temperature,  and  rj} 
is  the  stability  correction  terms  of  Panofsky  (Reference  42).  A  constant  value  of  0.1  m 
is  assumed  for  2„,  and  Zot  is  obtained  60m  a  relationship  ln{zofzot)  =  2  (Reference  41). 
Using  Equation  (A-24)  we  can  eliminate  T*  from  Elquation  (A-21)  to  obtain 
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R,  +  eRL  i  -taT%  +  m[©(20  -  Tg{PoI 


=  0 


(A  -  25) 


where 


m  =  kpaCpU^{l  ■¥  B  ^)Pr^[ln{{zi  + Zot)lzo} +2-xl)]  ^  (A  -  26) 


Po  is  a  reference  pressure  (1000  mb)  and  Pc  is  the  pressure  at  the  surface. 
Equation  (A-25)  may  be  linearized  by  noting  that 


(J2+1  _  j 

where  the  superscripts  n  and  n  +  1  denote  the  n  and  (n  +  1)‘*  time  siteps  of  integra¬ 
tion  (a  typical  time  increment  used  in  integration  is  1  minute).  After  substitution  into 
Equation  (A-25)  of  the  approximation 

(TS+')^  w4(TS)3T5+^-3(T5)^  ,  (A -27) 

we  obtain 


4e<7{TSy  +  miPo/PG)^^^^  + 


ISi 

Az, 


*  j 


pn+l 


=  (^)r;+>(i)  +  R,  +  efiu 

-f  3ecr(rQ)^  -h  m0’*(zi)  , 


(A  -  28) 


where  the  derivative  dTg/dzg  jo  is  replaced  by  a  forwEurd  finite-diflference  approximation 
( jn+i(i)_2’n+i)/^2,,T,(l)  is  the  soil  temperature  at  the  first  grid  level  from  the  surface, 
zmd  Aza  is  the  distance  between  the  surface  and  the  first  grid  level  in  the  soil  layer. 

Equation  (A- 16)  is  solved  numerically  in  finite-difference  form  by  Laasonen  (Ref¬ 
erence  11,  p.  189).  By  this  method  Equation  (A-16)  reduces  to  AT,  =  B  where  A  is  a 
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tridiagonal  matrix  and  B  is  a  column  vector.  The  solution  is  conveniently  obtained  by 
using  the  relation  (Reference  11,  p.  198) 


{Ts)t  =  Ei{T,)t+i+F^  ,  (A -29) 

where  {Ta)(  is  the  soil  temperature  at  the  grid  level  from  the  surface.  Expressions  for 
Et  and  Ft  when  i  >  1  eire  determined  from  the  finite-difference  form  of  Equation  (A- 16), 
and  Equation  (A-28)  determines  Ei  and  Fi.  From  Equations  (A-28)  and  (A-29),  we  obtain 


E^ 


-(&> 


4ea(TG)^  +  m 
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and 


(A  -  30) 


l-R,  +  fiii.  i  +3<(7(T^)'‘  +  me"(ii )]  f4_qi\ 

Numerical  integration  of  Equation  (A-16)  by  use  of  Equation  (A-29)  to  Equa¬ 
tion  (A-31)  is  rapid  since  no  iteration  is  required. 

The  incoming  direct  solar  radiation  flux  to  an  inclined  surface  is  obtained  from 
Kondratyev  (Reference  43); 


where 


Ra  =  Ro[A-\-  B  cos  ft  -t-  C  sin  ft] 


(A  -  32) 


A  =  cos  Q  sin  sin  ^  -I-  sin  a[cos  'J’„(tan  #  sin  $  sin  ^  —  sin  S  sec  $)]  (A  —  33a) 
B  =  cos  a  cos  $  cos  ^ -f  sin  a  cos  $  sin  cos  ^  ,  (A  —  336) 

and 

C  =  sin  a  cos  S  sin  #n  •  (A  —  33c) 

In  the  above  expressions  Ro  is  the  neeir  surface,  direct  solar  rsuliation  flux;  ft  is  the  solar 
hour  zmgle,  positive  clockwise  from  apparent  noon;  $  is  the  latitude;  S  is  the  declination 
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of  the  sun;  q  is  the  angle  of  incUnation  of  the  surface  relative  to  the  horizontal  plane;  and 
is  the  azimuth  of  the  projection  of  the  normal  to  the  surface  on  the  horizontal  plane, 
as  counted  from  the  plane  of  the  meridian  (azimuth  is  considered  positive  when  counted 
clockwise).  Since  the  maximum  change  in  the  solar  declination  6  in  24  hours  is  less  than  0.5 
degrees,  6  is  assumed  to  be  constant  during  a  given  day.  Spencer  (quoted  in  Reference  44, 
p.  63)  provides  a  formula  to  compute  8  in  radians, 


8  =  0.006918  -  0.399912  cos  +0.070257  sin 
-  0.006758  cos  2^0  +  0.000907  sin  2^0 

-0.002697  cos  3^0 +  0.001480  sin  3^0  ,  (34) 


where  the  angle  in  radians  is  related  to  the  Julian  day  by 


Bo 


27r(  Jrf  -  1) 

365 


{A  -  35) 


Equation  (A- 34)  estimates  8  with  a  maximum  error  of  0.0006  radians.  Solar  hour  f2  can  be 
obtained  if  the  longitude,  clocktime,  and  the  equation  of  time  are  known.  The  equation  of 
time  is  the  difference  between  the  local  apparent  time  and  a  fixed  mean  soIeu  time,  which 
is  derived  from  the  motion  of  a  celestial  equation  at  a  rate  equal  to  the  average  movement 
of  the  sun.  The  soIeu  hour  angle  Q  is  given  in  radians  by 


(A  -  36) 


where  is  the  true  solar  time  (local  apparent  time)  in  hoius.  The  true  solar  time  is 
obtained  from 


—  ic.t.  +  ^^long  +  tej  > 


(A  -  37) 


where  fc.t.i  ^^long,  and  are  the  clocktime,  the  longitude  correction,  and  the  equation 
of  time,  respectively.  The  longitude  correction  Eiccoimts  for  the  difference  between  the 
IoceJ  meridiEin  and  a  standard  meridian,  and  is  positive  if  the  local  meridian  is  esist  of  the 
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staxidard.  The  equation  of  time  is  provided  by  Spencer  (quoted  in  Reference  44,  p.  63)  as 
follows; 


t  =  —(0.000075 +  0.001868  cos  00 

TT 

-  0.032077  sin  0^  -  0.014615  cos  20^ 

-0.40849  sin  20o)  ,  (A -38) 

where  t^q  is  in  hours  and  0^  is  defined  by  Equation  (A-35).  Equation  (A-38)  has  a  maodmum 
error,  compared  with  values  tabulated  in  the  National  Almanac,  of  35  s  in  time. 

The  amoimt  of  solar  radiation  reaching  the  surfax:e  is  much  less  than  that  at  the 
top  of  the  atmosphere  due  to  many  factors,  including  molecular  scattering  and  absorption 
by  permanent  gases  such  as  oxygen,  ozone,  and  carbon  dioxide.  The  effect  is  parameterized 
by  Atwater  and  Brown  (Reference  45),  who  modified  the  original  form  by  Kondratyev 
(Reference  46)  to  include  the  effect  of  the  forward  Rayleigh. scattering.  The  expression  is 


G  =  0.485 
+  0.515 


L.041-0.16^- 


000949P  + 0.051 


cosZ 


■) 


{A  -  39) 


where  P  is  pressure  in  mb.  Other  important  factors  that  also  modify  the  amoimt  of  incom¬ 
ing  solar  radiation  include  water  vapor,  clouds,  and  airborne  peuticles.  Parameterizations 
for  these  factors  are  not  included  in  the  present  model. 

Currently  Ro  in  Equation  (A-32)  is  cedculated  fix>m 


i?o  —  RooG  ,  (A  —  40) 

where  Roo  is  the  incoming  radiation  flux  at  the  top  of  the  atmosphere  and  G  is  given  by 
Equation  (A-39). 

The  zenith  angle  Z  in  Equation  (A-39)  is  determined  fix>m  the  following  formula 
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cos  Z  =  sin  $  sin  b  +  cos  $  cos  ScosQ  .  (.4  —  41) 

Finally,  Ri  1,  the  long- wave  incoming  radiation  at  the  surface,  is  computed 
according  to  the  following  formula 

Rl  i=  Ro  i  cos  a  (A  -  42) 

where  Ro  1  is  the  long-wave  incoming  radiation  normal  to  horizontal  sin^ace  and  a  is  the 
angle  of  inclination  of  a  sloped  surface  given  by 

dz  dz 

Boimdary  conditions  for  U,V,Q,Q„,q,  and  £  along  the  upper  computational 
boimdary  are 

(U,V)  =  (l7„y,),  (A -44a,  6) 

where  Ug  and  Vg  are  geostrophic  wind  components  defined  as 

(C^.,v;)  =  (l//)(^),(f)  .  u-45a,6) 

Potential  temperature  and  the  mixing  ratio  are  specified,  and  tmrbulence  is  assumed  to 
vanish  along  the  upper  boundary.  Soil  temperature  at  30  cm  below  the  sirrface  is  also 
specified. 

The  lateral  boimdary  values  for  U,  V,  0,  Q„,  and  i  are  obtained  by  integrating 
the  corresponding  governing  Equations  (A-2),  (A-3),  (A-10),  (A-11),  (A-12),  and  (A-13), 
except  that  variations  in  the  horizontal  directions  are  all  neglected.  Variables  V,  and 
(^l  are  smoothed  at  each  time  step  by  using  the  values  at  four  neighboring  points,  i.e., 

=  (1  ~  +  0.25A($,+i,j  -I-  +  ^i,i+i)  ,  {A  -  46) 
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where  $  represents  either  i7,  V,g^,  or  q'^£,  and  A  =  0.5  is  used.  A  similar  expression,  but 
using  only  three  neighboring  points,  is  applied  to  the  values  at  the  lateral  boimdaries. 

Two  modifications  have  been  made  for  use  on  a  microcomputer.  First,  changes 
have  been  made  in  the  code  that  destroy  the  effects  of  terrain  slopes  at  the  boundaries. 
This  makes  it  possible  to  put  the  boimdaries  in  very  complex  terrain  without  constructing 
an  artificial  apron  around  the  area  of  interest.  The  second  change  was  to  force  the  first 
grid  vertical  cell  to  be  at  4  meters  abovegroimd.  This  change  was  to  ensure  that  the  lowest 
layer  was  sufliciently  near  the  groimd  to  generate  the  appropriate  slope  flows. 

3.  Initial  Values 

An  initial  wind  profile  at  the  southwestern  comer  of  the  computational  domain 
is  first  constmcted  by  assuming  a  logarithmic  variation  (initially  it*  =  0.2  m/s,  and  Zq  = 
0.1  meters)  from  the  ground  up  to  the  level  where  the  wind  speed  reaches  an  ambient  value 
(geostrophic  wind).  Initial  wind  profiles  at  other  grid  locations  are  obtained  by  scaling  the 
southwestern  comer  winds  to  satisfy  the  mass  continuity. 

The  vertical  profile  of  potential  temperature  is  initially  assumed  to  increase  lin¬ 
early  with  height.  Initial  potential  temperatures  are  Eissumed  to  be  uniform  in  the  horizon¬ 
tal  directions.  Initial  values  for  water  vapor  are  constructed  by  using  the  initial  potential 
temperature  profiles,  pressure  at  the  top  of  computational  domain,  and  observed  relative 
humidity.  The  turbulence  kinetic  energy  and  length  scale  are  initialized  by  using  the  ini- 
tizJ  wind  smd  temperature  profiles  and  the  relationships  resulted  from  the  Level  2  model. 
These  expressions  are  already  given  by  Yamada  (Reference  19)  and  are  not  repeated  here. 

4.  Numerical  Integrations 

The  partial  differential  Equations  (A-2),  (A-3),  and  (A-10)  to  (A-13)  are  inte¬ 
grated  by  using  the  Alternating  Direction  Implicit  (ADI)  method,  and  a  time  increment 
is  chosen  to  satisfy  Courant-Friedrich-Lewy  criteria.  To  increase  the  accuracy  of  finite- 
difference  approximations,  mean  and  turbulence  variables  are  defined  at  grids  that  are 
staggered  botn  in  horizontal  amd  vertical  directions.  Mean  winds,  temperatures,  and  wa- 
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ter  rapor  vary  greatly  with  height  near  the  surface.  In  order  to  resolve  these  variations, 
nonuniform  grid  spacings  axe  vised  in  the  vertical  direction, 

B.  RAPTAD  (RAndom  Particle  Transport  And  Diffusion.) 

A  brief  description  of  the  RAPTAD  model  is  given  here.  Locations  of  particles 
are  computed  from 


Xi{t  +  At)  =  Xi{t)  +  UpiAt  ,  {A  —  47) 

where 

Upi  =  Ui  +  Ui  ,  (A-  48) 

Q 

Ui(t  -f  At)  =  aui{t)  +  bouX  +  <5, -3(1  -  a)tLi.-  “  49) 

a  =  eip(-At/fx,*..)  ,  (A  -  50) 

and 

6  =  (l-a2)*/2  .  (A -51) 

In  the  above  expressions,  Upi  is  the  pzuticle  velocity  in  direction,  Ui  mean  velocity, 
Ui  turbulence  velocity,  ^  a  random  number  from  a  Gaussian  distribution  with  zero  mean 
and  unit  variance,  t£,j..  the  L£igrangian  integral  time  for  the  velocity  Uj,  cr„.  variance 
of  velocity  fluctuation  Uj,  and  ^,3  is  the  Dirac  delta.  The  last  term  on  the  right-hand 
side  of  Equation  (A-49)  was  introduced  by  Legg  and  Raupach  (Reference  47)  to  correct 
accumulation  of  particles  in  the  low  energy  areas.  The  mean  velocity  Ui  and  vertical 
velocity  variance  <t»,.  are  obtained  form  the  hydrodynamic  model  results.  The  Lagrangian 
time  scales,  =  20s,  =  5000s,  and  =  5000s,  are  used  in  this  study. 

In  the  previous  studies  (Reference  28;  Reference  39),  the  concentration  at  a  given 
time  and  location  was  determined  by  counting  the  number  of  particles  in  an  imaginary 
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sampling  volume.  The  computed  concentration  level  could  vary  considerably,  depending 
on  the  size  of  the  sampling  voliune  and  number  of  particles  used  in  the  computation.  For 
example,  if  the  sampling  volume  is  very  small,  the  concentration  distribution  becomes 
noisy.  On  the  other  hand,  if  the  sampling  volume  is  too  large,  the  concentration  distribu¬ 
tion  will  be  oversmoothed  (Reference  48).  Theoretically,  the  sampling  voliime  problem  is 
ehminated  by  releasing  an  infinite  number  of  particles  in  the  computation.  Of  course,  it  is 
impossible  in  practice,  or  at  least  very  expensive,  to  release  an  infinite  niimber  of  particles. 

A  “kernel”  density  estimator  is  used  in  this  study  where  each  particle  represents 
a  center  of  a  puff.  Various  functional  forms  can  be  assumed  to  express  the  concentration 
distribution  in  the  puff.  One  of  the  simplest  ways  is  to  assume  a  Gaussian  distribution 
where  variances  are  determined  as  the  time  integration  of  the  velocity  variances  encotm- 
tered  over  the  history  of  the  puff.  The  concentration  level  at  a  given  time  and  space  is 
determined  as  the  sinn  of  the  concentrations  each  puff  contributes.  The  kernel  method 
requires  no  imaginary  sampling  volximes  and  produces  smooth  concentration  distribution 
with  a  much  smaller  number  of  particles  than  required  for  the  previous  particle  method 
(Reference  47). 

Concentration  x  at  (X,  V,  Z)  is  estimated  by  using  the  following  expression: 


xiX,Y,Z) 


Q^t  ^  1  f  i(xfc-^\ 

(2;r)3/2  ^  2  al^  ) 


•  exp 


( 


Hyk-Yy\ 

2  <  ) 


1  {zk  -  Zf 

2  < 


l{z,  +  Z-2z,)^\ 

2  <  J 


(A  -  52) 


where  (x*,  y*,  z*)  is  the  location  of  k  th  particle;  and  are  standard  deviations  of 

a  Gaussian  distribution;  and  Zg  is  the  groimd  elevation.  The  variances  are  estimated  based 
on  Taylor’s  homogeneous  diffusion  theory  (Refer  ence  49).  For  example,  is  obtained  from 
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rl  =  2cl  f  f  mOHdt 

Jo  Jo 


=  2<7^iL,  t  +  tiyexp 


(-i)- 


<L 


y  /  ’ 


(A  -  53) 


where  a  correlation  function  R{C)  =  exp(  used.  Equation  (A-53)  is  approximated 

by 


and 


for  i  ) 

O’y  =  2tLyCrlt  for  t  >  2tiy  . 


{A  —  54a) 

(A  —  545) 


Although  the  turbulence  field  tinder  the  study  is  not  homogeneous,  we  assume  the  theoiy 
can  be  applicable  over  a  short  time  period,  such  as  an  integration  time  step  (10  seconds  in 
this  study).  Therefore, 


<T,(<  + At)  =  <Ty(t)  +  <7„A<  for<<2<i, 

and 

o-y(t  +  At)  =  <ry(t)  +  2tiyO-^At  for  t  >  2tx,y 

eire  used  in  this  study. 

In  a  similair  fashion, 


(A  -  55a) 

(A  -  556) 


<Tx(t  +  At)  =  <r*(t)  +  <7«At  for  t  <  2tii  , 
crl{i  +  At)  =  or*(t)  +  2tLx<7-^At  for  t  >  2tL*  , 
(Tj(t  +  At)  =  o',(t)  +  OttjAt  for  t  <  2t£,  , 


(A  —  56a) 
(A  -  566) 
(A  -  57a) 
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(rl{t At)  =  (xl{t) +  2tLza‘l,At  ioTt>2iLz  ,  {A  — 57b) 

where  the  standaxd  deviations  and  at  each  particle  location  axe  obtained  by 

interpolating  grid  values  of  a  computation  grid  volume  in  which  a  particle  is  located. 
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NUMERICAL  PROCEDURES 


A.  FINITE-DIFFERENCE  EQUATIONS 

The  prognostic  Equations  (A-2),  (A-3),  (A-10),  (A-11),  (A-12),  and  (A-13)  may 
be  expressed  in  the  following  general  form; 


where 


and 


4^{U,V,et,Qu„q\qH)  . 


{A  -  58) 

(A  -  59a) 
(A  -  59b) 
{A  -  59c) 

(A  -  60) 


The  coefficient  Ki  represents  the  horizontal  eddy  viscosity  coefficients  Kx  or  Kxy.  Simi¬ 
larly,  A’2  represents  Ky  or  Kxy,  and  K3  represents  the  vertical  eddy  viscosity  coefficients 
Km  or  Kh-  The  fourth  and  fifth  terms  on  the  right-hand  side  of  Equation  (A-58)  represent 
the  variable  to  which  the  equation  applies  and  external  forcing  functions,  respectively.  Ta¬ 
ble  A- 1  summarizes  Ki,K2,Kz,  A,  and  F for  the  prognostic  equations  for  U,  V,  0/,  Q w ,  5^ , 
and  q^£. 

The  ADI  method,  developed  by  Peaceman  and  Rachford  (Reference  50),  has 
second-order  accturacy  for  both  space  and  time  derivatives  and  is  unconditionally  stable. 
The  ADI  scheme  has  been  extensively  and  successfully  applied  in  the  simulation  of  various 
one-  and  two-dimensional  fluid  dynamics  problems  (Reference  51,  p.  95).  Generalization 
of  the  scheme  to  a  three-dimensional  space,  however,  reqmres  special  consideration,  as 
pointed  out  by  Richtmyer  and  Morton  (Reference  11,  p.  212);  otherwise,  the  xmconditiond 
stability  is  lost  and  accuracy  drops  to  0(At)  +  0[(Ax)^]. 
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TABLE  A-1.  COEFFICIENTS  ATj,  ATa,  A,  AND  F  In  EQUATION  (A-58) 


<!> 

Ki 

K2 

Kz 

A 

U 

Kz 

Kry 

Km 

0 

V 

Kzy 

Ky 

Km 

0 

A0/ 

K, 

Kh 

Kh 

0 

Qw 

Ky 

Kh 

0 

al 

2 

K, 

Ky 

qe'St 

9/A1 

qH 

Kz 

Ky 

qi^t 

qj  Ai  X 

1 + (k 


_ F _ 

H  ( 1 
H-2,  I  pCp  ar*  ^  dz*  J 

0 

+  ^1^)  + 

^Fi  -  ^0)  +  ^9w0r 


The  finite-difference  version  of  Equation  (A-58)  may  be  written  according  to  ADI 
method  as 


=  iA.(r  +  D  +  A,^»  +  A,r  -Ar+F  ,  (A-  61) 

+  r)  +  +  «>”)  +  -Ar  +  F  ,  (A-  82) 

and 

-6^  1  1 
~t  = 

+  +  ^“)  -Ar  +  F  ,  (A  -  63) 

where  <^"  and  axe  values  at  t  =  nAt  emd  (n  -f  l)Af,  respectively,  and  At  is  a 

time  increment  of  integration.  The  value  represented  by  <f>*  and  <f>**  are  the  first  and 
second  intermediate  values  necessary  to  obtEun  firom  Richtmyer  and  Morton 
(Reference  11,  p.  213)  pointed  out  that  if  <l>*  in  the  vertical  derivative  term  of  Equation  (A- 
63)  is  replaced  by  the  most  recent  value  then  the  \mconditional  stability  is  lost.  The 
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Eimplification  factor  for  the  complete  cycle  given  by  Equations  (A-61)  to  (A-63)  is  shown 
to  be  not  greater  than  1,  indicating  unconditional  stability  (Reference  11,  p.  213). 

Equations  (A-62)  and  (A-63)  may  be  simplified  by  subtracting  Equation  (A-61) 
from  Equation  (A-62)  and  subtracting  Equation  (A-62)  from  Equation  (A-63),  respectively. 
The  resulting  equations  are 


and 


A**  -  A*  1 

- - —  =  -ArU**  -  <^”) 

At  2  ^  ’ 


-  6**  1  , , 

At  2  ^  ’ 


Equations  (A-61),  (A-64),  and  (A-65)  may  be  rewritten  as 


{A  -  64) 


(A  -  65) 


(l  +  AAt  -  YA,)r  =  +  F}  (yl  -  66) 

(l- YA.<1”  , 


^"+'  =  4,-  -  YA  . 


(A  -  67) 


(A  -  68) 


Equations  (A-66)  and  (A-68)  may  be  further  transformed  to  the  general  form 


-h  —  C(<f>i+i  =  Df  ,  (-4  “  69) 

where  i  =  i,j,k  represents  the  i,j,  and  k  th  grid  points  in  the  x,y,  and  z  directions, 
respectively.  The  space  derivatives  A*^,  Ay^,  and  Ag<f>  are  approximated  by  centered 
finite  differences.  Coefficients  At,  Bt,  Ct,  and  Dt  may  be  identified  by  comparing  term  by 
term  the  expanded  forms  of  Equations  (A-66)  to  (A-68)  with  those  of  Equation  (A-69). 
The  results  are  given  in  Table  A-2.  The  coefficients  ai,bj,  and  cj  are  dXfdx,dYfdy,  and 
dZjdz* ,  respectively,  where  {X,Y,Z,)  are  grid  coordinates,  whose  increments  are  one. 
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TABLE  A-2.  COEFFICIENTS  Ae,Bt,Ct,  and  Di  FOR  EQUATION  (A-69). 


Ak 

(cibA(/2)| 

^  +  (ciir3)t_i/2| 

Ck 

(ciAr/2) 

1  — 2^  +  (cA'3)it-H/2| 

Bk 

1  -f  Ak  -|-  Cfc  -t-  AAt 

Dk 

-1-  A^<^”  Aj,<^"  -H  F) 

A, 

(a,At/2)| 

^A{aKi)i-\l2^ 

c. 

(a.Ai/2)| 

^+iaKi)i  +  l/2\ 

Bi 

1  +  A,-  -|-  < 

<  } 

’•'i 

Di 

<t>*  -  (At/2)A,<ji" 

(6yA(/2)| 

-f  (6Ar2)_,_i/2| 

Cj 

(6yA(/2)| 

I  +  (^■^2)j-(-1/2| 

1  -f-  A  j  -f  Cj 

Di 

<i>y  -  {At/2)Ayr 

Equation  (A-69)  may  be  expressed  in  the  general  form  where  ^  is  a 

tridiagonad  matrix  whose  elements  are  given  by  At,  Bt,  and  Ce,  and  ^  is  a  column  vector 
whose  elements  are  given  by  Dt.  Solutions  for  Equation  (A-69),  whose  coefficient  matrix 
is  tridiagonal,  are  obtained  by  a  direct  elimination  method  (Reference  11,  p.  200). 

Accuraw:y  of  a  finite-difference  approximation  is  enhanced  if  more  grid  points 
are  placed  where  the  variables  vary  greatly  with  the  co<x:dinates.  For  example,  wind, 
temperature,  and  water  vapor  are  known  to  vary  considerably  with  height  in  the  surf2M:e 
layer.  Those  variables,  however,  vary  much  more  slowly  with  height  in  the  layers  away 
&om  the  surface.  Thus,  vertical  grids  are  spared  linearly  near  the  surfsM:e  and  parabolically 
away  from  the  surface: 
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z  =  c\k  for  k  <  kt 


{A  -  70a) 


z  =  Cik  +  C2k  +  C3  for  k  >  kt  ~  706) 

where  k  is  the  vertical  grid  level  whose  increment  is  one  and  kt  is  a  transition  level  between 
the  linear  and  parabolic  variations.  Constraints  ci,C2,C3  and  C4  are  determined,  depending 
on  the  problem  to  be  solved. 

In  order  to  increase  the  accuracy  of  finite-difference  approximations  and  to  sup¬ 
press  computational  noise  whose  wave  length  is  2A,  where  A  is  a  grid  increment,  mean 
2ind  turbulence  variables  are  defined  on  a  staggered  grid. 

Figure  A-1  shows  relative  locations  of  variables  in  a  grid  volume:  U  and  V  in 
Figure  A-la;  and  Quj  in  Figure  A-lb;  and  and  in  Figure  A-lc. 

B.  BOUNDARY  CONDITIONS 

As  shown  by  Mellor  (Reference  16),  oxir  model  simulated  the  surface  layer  data 
reported  by  Businger  et  ed.  (Reference  52)  quite  well.  A  rather  coarse  vertical  grid  spacing 
is  necessary  for  three-dimensional  modeling  due  to  limitations  of  storage  and  time.  Thus, 
in  order  to  increase  accuracy,  empirical  surface-layer  paurameterization  is  used  for  the  fii-st 
and  second  grid  levels  above  the  ground.  Wind,  temperature,  and  water  vapor  profiles  in 
the  surface  layer  may  be  approximated  by  the  following  formulas: 

*  =  t[^«{(^  +  ^o)Ao}  -  V’m(C)]  ,  {A  -  71) 

Ua  K 

^  ^  _ 72) 

g.L^bq»(P)  =  _  73) 

,^(r)  =  B?/\2(^„-<)V3  .{74) 

and 
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q^i{z)  =  kzq^{z)  . 


{A  -  75) 


In  the  above  equations,  ^(2),0£(2),  and  Q^i^)  are  the  abbreviations  for  2,t), 

Qt{x,y,z,t),  and  Quj{x,y,z,t),  respectively,  and  |  |  is  the  horizontal  wind  speed  ({7^  + 

y2)i/2  'Perms  u*,T,,  and  Q.,  defined  as  u*  =  y/rjp,  T,  =  H f pCpU^^  and  Q»  =  Ej pu*, 
are  friction  velocity  and  scales  for  temperature  and  water  vapor,  respectively.  Terms  p,  H, 
and  E  are  surface  stress,  total  sensible  heat,  and  rate  of  evaporation,  respectively;  total 
sensible  heat  being  defined  as  H  =  wdy  =  (1  +  O.61Qu,)id0  +  0.61uJ^.  The  parameter 
Q  is  a  nondimensional  height  z/L,  where  L  is  the  Obukhov  length  -  ul/k^gH;k  is  the 
von  Karman  constant;  Zo,Zot,  and  Zov  are  roughness  lengths  for  wind,  temperature,  and 
water  vapor,  respectively;  and  Pr  and  Sc  are  the  turbtilence  Prandtl  and  Schmidt  numbers, 
respectively.  Terms  V’m,  V’A?  ^nd  (Reference  42)  are  correction  terms  for  the  atmospheric 
stability,  and  their  functional  forms  are  given  by 


(A  -  76) 

(A  -  77) 

and 

1 

II 

(A  -  78) 

where  <l>m,<f>h,<^v  are  nondimensional  wind,  temperatiure,  and  water  vapor  gradients,  re- 

spectively.  The  following  formulations  (Reference  53;  Reference  10) 
emd  <j>v  imder  unstable  conditions, 

are  used  for  <i>m-,4>hi 

(A  -  79) 

(A  -  80) 

and 

(A  -  81) 
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For  stable  conditions, 


MO  =  MO  =  MO  =  1  +  5C  ■  {A-  82a,  6,  c) 

The  roughness  lengths  Zo,  Zot,  and  Zov  are  specified  over  land  and  are  computed  over  water 
according  to  the  following  formulas; 


and 


a 


{A  -  83a) 
(A  -  836) 


Zov  =  -p~~  (A  —  83c) 

tcxt^ 

(see  Reference  54;  Reference  55),  where  7  is  a  constant  whose  numerical  value  varies  from 
0.008  to  0.032  for  water  surfaces  vzuying  from  a  lake  to  a  rough  ocean  and  a  and  h  axe  the 
molecular  diffusivities  for  heat  and  water  vapor,  respectively.  Integrations  of  Equations  (A- 
76)  to  (A-78)  may  be  performed  easily  (Reference  55).  The  scales  u,,r,,  and  Q*  may  be 
obtained  by  solving  iteratively  Equations  (A-71)  to  (A-81)  omder  rmstable  conditions. 
Solutions  may  be  obtained  without  iteration  for  stable  conditions  over  land.  Iterations  axe 
required  over  the  water,  however,  since  the  roughness  heights  axe  also  functions  of  u*,  eis 
seen  in  Equations  (A-82a-c.) 

Solutions  of  Equation  (A-69)  may  be  obtained  using  the  bovmdaxy  conditions 
discussed  above  and  following  the  algorithms  described  by  Richtmyer  and  Morton  (Refer¬ 
ence  11,  p.  200).  Solutions  axe  eissumed  to  be  given  by  a  relation 


<l>t  =  +  Et 

where  Ee  and  Ft  axe  determined  from 


(A  -  84) 


Et 


Ce 

Bt  —  AfEt-i 


,e>i 


(A  -  85a) 
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_  -Df  +  AtFt-i 
'  Bi-AiEi-i  ’  ^  ’ 


(A  -  85b) 


where  coefficients  Ai,Bt,Ct,  and  De  axe  given  in  Table  A-2.  The  first  values  Ei  and  Fi 
may  be  determined  from  the  surface  boimdary  conditions  given  by  Equations  (A-71)  and 
(A-75).  For  example,  for  Y  ©<5 


tn{{Z2  +  Zo)lZo]  -  IpmiCz) 


(A  -  86) 


in{(zi  +  Zot)/zot}  -  rphiCi) 
£n{(z2  +  Zot)lzot}  -  tkhiCz) 


0/(0) 


(A  -  87) 


where  zi  and  Z2  are  the  first  and  second  grid  levels  above  the  surface,  respectively.  Com¬ 
parison  between  Equation  (A-84)  and  Equation  (A-86),  where  <j>  =  Xi  yields 


Ei  = 


tn{{zi  +^o)/go}-V>m(Cl) 

in{{z2  +  Zot)IZot}  -  i’hiCi) 


{A  -  88a) 


Fi  =  0  . 


(A  -  886) 


Similarly,  comparison  between  Equation  (A-84)  and  Equation  (A-87),  where  d  =  R  ,  yields 


El  = 


£n{{zi  -I-  Zot)/zot}  -  ^ft(Ci) 

£n{{z2  +  Zot)/zot}  -  ^hiCz) 


(A  -  89a) 


Fi=(l  -Fi)0,(O)  . 


(A  -  896) 


Similar  expressions  as  in  Equations  (A-89a)  and  (A-89b)  are  valid  for  Q^,  except  that 
Zot  is  replaced  by  Zov  When  boimdary  values  are  specified,  as  in  the  case  of  and 
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[Equations  (A-74)  and  (A-75)],  Ei  and  Fi  are  easily  obtained.  By  comparing  Equation  (A- 
S4)  with  Equation  (A-74)  we  obtain  for 


Ei=0  (A  -  90a) 

and 

Fi  =  .  (A  -  906) 

Similarly,  by  comparing  Equation  (A-84)  with  Equation  (A-75)  we  obtain  for 

El  =  0  (A  -  91a) 

and 

Fi=kziq^{zi)  .  (A  — 916) 

Once  El  and  Fi  axe  determined  as  described  above,  E/  and  Fi{£  >  1)  may 
be  computed  according  to  Equations  (A-85a)  and  (A-85b),  respectively.  Then  solutions 
may  be  computed  from  Equation  (A-84),  provided  the  value  at  the  upper  boimdary  is 
determined  from  upper  bovmdary  conditions  as  follows.  The  boundary  conditions  at  the 
top  are  assumed  to  be 


U  =  Ug(Hy,  V  =  Vg(H)  ,  (A -92a,  6) 

e  =  e(H)  ,  (A -93) 

Q^  =  QUH)  =  0  ,  (A -94) 

9^  =  0  ,  (A  -  95) 

and 

=  0  .  (A -96) 


At  lateral  boundaries,  wind,  temperature,  water  vapor,  g^,  and  g^£  are  computed 
from  the  one-dimensional  versions  of  the  corresponding  Equations  (A-2),  (A-3),  (A-10), 
(A-11),  (A-12),  and  (A-13). 
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